Extra parameter using iterative method in ODE45

I am not able to find out why the extra parameter/alpha is not converging to the exact solution. It just takes the first value in the interval.

 %
zspan=[0,400];
v0mat = [1 0.1 1];
tol = 0.0000001; % Treshold
MAX = 500; % Max Iteration
v2 = 20; % Initial value 
interval = [0.1 0.2]; % alpha interval ( alpha is supposed to be 0.116)
a = interval(1);
b = interval(2);
alpha = (a+b)/2;
zsol = {};
v1sol = {};
v2sol = {};
v3sol = {};
for k=1:size(v0mat,1)
  v0=v0mat(k,:);
  [z,v]=ode45(@(z,v)rhs(z,v,alpha),zspan,v0);
  zsol{k}=z;
  v1sol{k}=v(:,1);
  v2sol{k}=v(:,2);
  v3sol{k}=v(:,3);
  mssflx{k} =v(:,1);
  vel{k}=v(:,2)./mssflx{k};
end
iter = 1;
while(iter<MAX)
 alpha= (a + b)/2;
[z,v]=ode45(@(z,v)rhs(z,v,alpha),interval,v0);
for k12 = 1:length(vel)
  zsol04(k12) = interp1(real(vel{k12}), real(zsol{k12}),0.4);
end
if abs(zsol04-v2) > tol
 b=alpha;
 else
 a = alpha;
 end
if abs(zsol04-v2)< tol
     disp('Sucess');
          break  
 end
iter = iter + 1;
end
for r=1:length(v2sol)
q(r)=r;
end
figure
scatter(q,zsol04,'g')
function parameters=rhs(z,v,alpha)
            db= 2*alpha-(v(1).*v(3))./(2*v(2).^2);
            dw= (v(3)./v(2))-(2*alpha*v(2)./v(1));
            dgmark= -(2*alpha*v(3)./v(1)); 
   parameters=[db;dw;dgmark];
    end

Answer 1

[Walter Roberson]

You have

if abs(zsol04-v2) > tol
 b=alpha;
 else
 a = alpha;
 end

That logic is wrong. You are obviously trying code bisection, so after you have checked that you are not within tolerance, you should be testing the sign of zsol04-v2 to determine whether to move a or move b .

[Dereje]

Like you said I have corrected it:

 while abs(zsol04-v2) > tol
 alpha= (a + b)/2;
 % After some lines 
if zsol04-v2 > 0
 b=alpha;
 else
 a = alpha;
end
iter = iter + 1;
if (iter > MAX)
     return;
  end

But still no change.

[Walter Roberson]

You have

[z,v]=ode45(@(z,v)rhs(z,v,alpha),interval,v0);
for k12 = 1:length(vel)
  zsol04(k12) = interp1(real(vel{k12}), real(zsol{k12}),0.4);
end

You are interpolating based upon vel and zsol, but neither of those are being changed by the ode45.

Question: what is the purpose of your z0mat ? You loop over the rows of it, but there is only one row. Your code

zsol04(k12) = interp1(real(vel{k12}), real(zsol{k12}),0.4);

is looping over the rows, but right after that loop you test all of zsol04,

if abs(zsol04-v2) > tol

which is equivalent of

if all(abs(zsol04-v2) > tol)

so it is as-if you are required to find simultaneous solutions. But then you

if zsol04-v2 < 0

which is equivalent to

if all(zsol04-v2 < 0)

and that is not the same as

if ~all(zsol04-v2 >= 0)

because some entries might be true and some might be false.

It might not be possible to find a single alpha value that satisfies all of the rows of v0mat simultaneously.

[Dereje]

The purpose of v0mat is to create a cell array(I was using many initial conditions in my previous case and took it as it is. It is right no use in this case ). Like you said the key is for the ode45 to change while vel and zsol changes and my problem was also that.