Why does my code for shooting method using 'ODE45' or 'ODE23s' does not converge to the boundary value.?

Question

naygarp 2018-5-1

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/398470-why-does-my-code-for-shooting-method-using-ode45-or-ode23s-does-not-converge-to-the-boundary-val

评论： Torsten 2018-5-4

采纳的回答： Torsten

在 MATLAB Online 中打开

I have used shooting method with ' ode45' or ' ode23s'.

But , the solution doesn't converge and it takes a lot of time.

The equations are

f"=g(g^2+gamma^2)/(g^2+lambda*gamma^2) ------ (1)
g'= (1/3)*f'^2-(2/3)*(f*f")+ Mn*f' ------------------------(2)
t"+Rd*t"+ 2*Pr*f*t'/3+ Nb*t'*p'+Nt*(t')^2= 0------(3)
p"+(2*Lew*f*p')/3+ Nt*t"/Nb= 0 ------------------------(4)

With the initial and boundary conditions

f(0)=0, f'(0)=1, t(0)=1, p(0)=1
f'(infinity)=0, t(infinity)=0, p(infinity)=0

The code for shooting method using ode45 is

function shooting_method
clc;clf;clear;
global lambda gama Pr Rd Lew Nb Nt Mn
gama=1;  
Mn=1;
Rd=0.1;
Pr=10;
Nb=0.2;
Lew=10;
Nt=0.2;
lambda=0.5;
x=[1 1 1];
options= optimset('Display','iter');
x1= fsolve(@solver,x);
end
function F=solver(x)
options= odeset('RelTol',1e-8,'AbsTol',[1e-8 1e-8 1e-8 1e-8 1e-8 1e-8 1e-8]);
[t,u]= ode45(@equation,[0 10],[0 1 x(1) 1 x(2) 1 x(3)],options)
s=length(t);
F= [u(s,2),u(s,4),u(s,6)];
%deval(0,u)
plot(t,u(:,2),t,u(:,4),t,u(:,6));
end
function dy=equation(t,y)
global lambda gama Pr Rd Lew Nb Nt Mn
dy= zeros(7,1);
dy(1)= y(2);
dy(2)= y(3)*(y(3)^2+gama^2)/(y(3)^2+lambda*gama^2);
dy(3)= y(2)^2/3-(2*y(1)*y(3)*(y(3)^2+gama^2))/(3*(y(3)^2+lambda*gama^2))+Mn*y(2);
dy(4)= y(5);
dy(5)= -(2*Pr*y(1)*y(5))/(3*(1+Rd)) - (Nb*y(5)*y(7))/(1+Rd) - (Nt*y(5)^2)/(1+Rd);
dy(6)= y(7);
dy(7)= -((2*Lew*y(1)*y(7))/3)+ (Nt/Nb)*((2*Pr*y(1)*y(5))/(3*(1+Rd)) + (Nb*y(5)*y(7))/(1+Rd) + (Nt*y(5)^2)/(1+Rd));
end

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Torsten 2018-5-2

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/398470-why-does-my-code-for-shooting-method-using-ode45-or-ode23s-does-not-converge-to-the-boundary-val#answer_318214

Try to start with the solution you get from "bvp4c" for the vector x.

Best wishes

Torsten.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 2

Jan 2018-5-2

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/398470-why-does-my-code-for-shooting-method-using-ode45-or-ode23s-does-not-converge-to-the-boundary-val#answer_318233

编辑：Jan 2018-5-2

在 MATLAB Online 中打开

You want to get:

f'(infinity)=0, t(infinity)=0, p(infinity)=0

but you integrate on the interval [0, 10]. 10 is not infinity. It is possible, that there is no possible start value, which reaches the wanted final point at the time 10.

Another problem can be the standard limitation of the single shooting methods: if a certain parameter causes a trajectory with Inf or NaN values, convergence is impossible. Then a multiple-shooting approach can help. Ask Wikipedia for details.

You use ode45 or ode23s? One is for non-stiff, the other for stiff systems. Using them by trial and error seems to be a very relaxed method of applied mathematics.

8 个评论
显示 6更早的评论隐藏 6更早的评论

naygarp 2018-5-3

在 MATLAB Online 中打开

 function shooting_method
clc;clf;clear;
global lambda gama Pr Rd Lew Nb Nt Mn
gama=1;  
Mn=1;
Rd=0.1;
Pr=10;
Nb=0.2;
Lew=10;
Nt=0.2;
lambda=0.5;
x=[1 1 1];
options= optimset('Display','iter');
x1= fsolve(@solver,x);
end
function F=solver(x)
options= odeset('RelTol',1e-8,'AbsTol',[1e-8 1e-8 1e-8 1e-8 1e-8 1e-8 1e-8]);
sol = ode15s(@equation,[0 4],[0 1 x(1) 1 x(2) 1 x(3)],options)
%sol= [t,u];
s=length(t);
F= [u(s,2),u(s,4),u(s,6)];
y1 = deval(u(0,:))
plot(t,u(:,5),t,u(:,7));
end
function dy=equation(t,y)
global lambda gama Pr Rd Lew Nb Nt Mn
dy= zeros(7,1);
dy(1)= y(2);
dy(2)= y(3)*(y(3)^2+gama^2)/(y(3)^2+lambda*gama^2);
dy(3)= y(2)^2/3-(2*y(1)*y(3)*(y(3)^2+gama^2))/(3*(y(3)^2+lambda*gama^2))+Mn*y(2);
dy(4)= y(5);
dy(5)= -(2*Pr*y(1)*y(5))/(3*(1+Rd)) - (Nb*y(5)*y(7))/(1+Rd) - (Nt*y(5)^2)/(1+Rd);
dy(6)= y(7);
dy(7)= -((2*Lew*y(1)*y(7))/3)+ (Nt/Nb)*((2*Pr*y(1)*y(5))/(3*(1+Rd)) + (Nb*y(5)*y(7))/(1+Rd) + (Nt*y(5)^2)/(1+Rd));
end

Error message that I receive is

    Undefined function or variable 't'.
Error in shooting_method_scientific_rana_ode45>solver
(line 21)
s=length(t);
Error in fsolve (line 230)
            fuser = feval(funfcn{3},x,varargin{:});
Error in shooting_method_scientific_rana_ode45 (line
14)
x1= fsolve(@solver,x);
Caused by:
    Failure in initial objective function evaluation.
    FSOLVE cannot continue.

Torsten 2018-5-4

在 MATLAB Online 中打开

function F=solver(x)
options= odeset('RelTol',1e-8,'AbsTol',[1e-8 1e-8 1e-8 1e-8 1e-8 1e-8 1e-8]);
[t,u] = ode15s(@equation,[0 4],[0 1 x(1) 1 x(2) 1 x(3)],options)
F= [u(end,2),u(end,4),u(end,6)];
y1 = u(1,:);  % should be equal to [0 1 x(1) 1 x(2) 1 x(3)]
plot(t,u(:,5),t,u(:,7));
end