Extra parameter/alpha in ODE45

1 回表示 (過去 30 日間)
Dereje
Dereje 2018 年 4 月 25 日
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

回答 (0 件)

カテゴリ

Help Center および File ExchangeHolidays / Seasons についてさらに検索

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by