Extra parameter using iterative method in ODE45

1 次查看（过去 30 天）
Dereje
Dereje 2018-4-15
评论： Dereje 2018-4-16
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

请先登录，再回答此问题。

回答（1 个）

Walter Roberson
Walter Roberson 2018-4-15
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 .
  3 个评论
显示 1更早的评论
Walter Roberson
Walter Roberson 2018-4-15
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
Dereje 2018-4-16
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Colormaps 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by