How do I use a for loop in my ode15s based code for shooting method and get multiple graphs?

2 views (last 30 days)
I am using ode15s solver to solve a set of odes by shooting technique and obtain graphs of the solutions.
while trying to vary some parameters in the equations I have used a for loop.
But I am not getting the proper graphs.
How do I use the 'for' loop properly to get the accurate graphs.
Here in the code below I need to vary for three values of the 'lambda' parameter, so that I could obtain 3 different graphs in one figure
function shooting_method
clc;clf;clear;
global lambda gama Pr Rd Lew Nb Nt Mn
gama=1;
Mn=6;
Rd=0.1;
Pr=10;
Nb=0.3;
Lew=10;
Nt=0.3;
pp = [0.5 1 1.5];
for i=1:numel(pp)
lambda = pp(i);
%lambda=0.5;
x=[1 1 1];
format long
options= optimset('Display','iter');
x1= fsolve(@solver,x);
end
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] = ode15s(@equation,linspace(0,2),[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(:,2),'LineWidth',2);hold on %t,u(:,4));
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

Answers (1)

Torsten
Torsten on 29 May 2018
Edited: Torsten on 30 May 2018
function shooting_method
clc;clf;clear;
global lambda gama Pr Rd Lew Nb Nt Mn
gama=1;
Mn=6;
Rd=0.1;
Pr=10;
Nb=0.3;
Lew=10;
Nt=0.3;
pp = [0.5 1 1.5];
for i=1:numel(pp)
lambda = pp(i);
%lambda=0.5;
x=[1 1 1];
format long
options= optimset('Display','iter');
x1= fsolve(@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,linspace(0,2,20),[0 1 x1(1) 1 x1(2) 1 x1(3)],options)
U(i,:,:)=u;
T(i,:)=t;
end
plot(T(1,:),U(1,:,2),T(2,:),U(2,:,2),T(3,:),U(3,:,2))
end
  3 Comments
Torsten
Torsten on 1 Jun 2018
Of course,you will have to run the above modified function "shooting_method" together with your two functions "solver" and "equation" from above.

Sign in to comment.

Categories

Find more on Numerical Integration and Differential Equations in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!