Extra parameter/alpha in ODE45

1 次查看（过去 30 天）

Dereje 2018-4-25

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/397284-extra-parameter-alpha-in-ode45

I am not able to find out why the extra parameter/alpha is not converging to the exact solution. My code seems ok. I need your help.

    %
zspan=[0,400];
v0mat = [1 0.1 1];
tol = 0.0000001; % Treshold
MAX = 500; % Max Iteration
v2 = 10; % 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;
for k12 = 1:length(vel)
  zsol04(k12) = interp1(real(vel{k12}), real(zsol{k12}),0.4);
end
while abs(zsol04-v2) > tol
 alpha= (a + b)/2;
[z,v]=ode45(@(z,v)rhs(z,v,alpha),interval,v0);
if zsol04-v2 > 0
 b=alpha;
 else
 a = alpha;
end
iter = iter + 1;
if (iter > MAX)
     return;
end
end
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