Parameter calculation by using loops in ODE45

Dereje

2020 1 16

0 个回答

2 次查看（30 天）

0 个投票

For any initial values(Ofcourse with limits) I used it either excutes on the 4th iteration or it would be beyond the MAX iteration.

I am afraid something is wrong and I need a help.

zspan=[0,400];
v0mat = [1 0.05 1];
tol = 10^-4; % Treshold
MAX = 1000;
v2 = 0.4;
interval = [0.01 0.2]; % alpha interval 
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);
end
iter = 1;
while(abs((v(:, 2)) - v2) > tol)
 alpha= (a + b)/2;
[z,v]=ode45(@(z,v)rhs(z,v,alpha),interval,v0);
 if(abs(v(:,2))-v2 > 0)
b = alpha;
 else
 a = alpha;
 end
 
 iter = iter + 1;
 if(iter > MAX)
 return;
 end
end
for k=1:size(v0mat,1)
    figure(1)
     plot(v2sol{k},zsol{k},'g')
      xlabel('Velocity,w')
       ylabel('Height, z')
      grid on       
end
function parameters=rhs(z,v,alpha)
            Nsqr = 0.001;
            db= 2*alpha-(v(1).*v(3))./(2*v(2).^2);
            dw= (v(3)./v(2))-(2*alpha*v(2)./v(1));
            dgmark= -Nsqr-(2*alpha*v(3)./v(1)); 
   parameters=[db;dw;dgmark];
end

16 个评论
显示 14更早的评论隐藏 14更早的评论

Walter Roberson 2020-1-17

编辑：Walter Roberson 2020-1-17

In your context, timepoint is z value. Your ode45 call will process over z in zspan, producing one row for each of a number of z values. (abs(v(:,2))-v2 > 0) is a vector result that might be true for some of the z entries but false for other times. You need to decide which z the condition must hold for.