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

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

Answer 1

[Torsten]

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

Best wishes

Torsten.

Answer 2

[Jan]

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.

[Torsten]

Using them by trial and error seems to be a very relaxed method of applied mathematics.

But isn't "applied mathematics" the same as "trial and error" :-) ?

[naygarp]

Speaking of hit and trial, if I minimize the interval to [0 4] or [0 3] then convergence is obtained within a short time. But in bvp4c solver, taking a interval value of [0 10] or [0 20] etc. the convergence is arrived regardless.

[Torsten]

Maybe integration with ODE45 becomes more and more inexact over longer distances.

[Jan]

Maybe setting the options of fsolve helps to improve the convergence, see https://www.mathworks.com/help/optim/ug/fsolve.html#buta__s-options.

[naygarp]

One more thing, how do I use 'deval' to extract the values of all the 'u' at t=0. I tried several ways of using deval but everytime I get some error

[Jan]

Please post your code and a copy of the error message. Then it is possible to suggest an improvement.

[naygarp]

 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]

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