How do I use a for loop in my ode15s based code for shooting method and get multiple graphs?
조회 수: 1 (최근 30일)
이전 댓글 표시
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
댓글 수: 1
MINATI
2019년 2월 2일
Dear naygarp
Can you share the question or the paper of whose this is the code
minatipatra456@gmail.com
답변 (1개)
Torsten
2018년 5월 29일
편집: Torsten
2018년 5월 30일
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
Torsten
2018년 6월 1일
Of course,you will have to run the above modified function "shooting_method" together with your two functions "solver" and "equation" from above.
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!