using shooting method for coupled ode

Hi all,
please help,I'm looking to solve the following system of equations with boundary conditions using the shooting method:
The equation and m-file attached .
I've found the solution using the BVP4C solver but need to also be able to find the solution using the shooting method.
I'm really quite new to MATLAB and don't really know where to start!
Any help anyone can give me would be greatly appreciated. Thanks!

1 件のコメント

darova
darova 2020 年 4 月 13 日
Can you rewrite this part?
function f= projfun(x,y)
f= [y(2)
y(1)*y(3)+x*y(2)*y(3)+y(1)*y(2)+(y(9)*(1/alfa*Pr*x))-(y(2)/x)-(y(1)/(x*x))
(-1/x)*(4*y(3)+y(2)+(y(1)/x))
4*y(3)*y(3)+((x*y(3)+y(1)-(1/x))*((-1/x)*(4*y(3)+y(2)+(y(1)/x))))
y(6)
Pr*x*y(3)*y(6)+Pr*y(1)*y(6)-(y(6)/x)
y(8)
(1/B*Tinf)*(alfa*Pr*y(1)*y(8)+alfa*Pr*x*y(3)*y(8)-(((DT*Bs*Fiinf)/(Tinf*Dfi))*((Pr*x*y(3)*y(6)+Pr*y(1)*y(6)-(y(6)/x))-(y(6)/x))))
y(9)
];
end
To make it look like:
function f= projfun(x,y)
E
E = y(1);
dE = y(2);
F = y(3);
dF = y(4);
% ...
f = zeros(9,1);
f(1) = dF;
f(2) = long_expr1;
f(3) = long_expr2;
% ...
end
  1. Use variable instead of y(1),y(2) for long code
  2. Why not define f function as you did before?
example
[1 - 2; 2] % work
[1 -2; 2] % doesn't work
I can't check your function because it's difficult to read. I don't want to write new one

サインインしてコメントする。

回答 (1 件)

darova
darova 2020 年 4 月 14 日
編集済み: darova 2020 年 4 月 14 日

1 投票

Thank you more readable code
Here are some mistakes
Some note:
Shouldn't last equations be like this? What is g? There is no g in your equations
% g = y(7);
E = y(7);
dE = y(8);
% dy(7) = dg;
dy(7) = -1*(4*E+df+f); % E'
dy(8) = (4*E*E)+((E*x+f-1)*(-1*(4*E+df+f))); % E''
I also suggest you not to use global: LINK
You can just use nested functions
function proj
% global variables - no need
% code
function dy = projfun(x,y)
% global variables - no need
% code
end
end

14 件のコメント

T K
T K 2020 年 4 月 14 日
編集済み: darova 2020 年 4 月 14 日
OK,Doctor darova thanks and appreciation .
I have modified all the errors that you mentioned ,
It remains an error only for the Variable g', please help me .
How to write the definition of g' in the function f'' ???
E = y(7);
dE = y(8);
g = y(9); %is it coorect for the definition of g'??
dy(1) = df;
dy(2) = (E*f)+(E*x*df)+(f*df)+((alfa*Pr*dg/(rho)))-(df)-f; %is it coorect for definition of g'??
dy(3) = dt;
dy(4) = (Pr*dt*E)+(Pr*f*dt)-(dt);
dy(5) = dc;
dy(6) = (1/(B*Tinf))*((alfa*Pr*dc*f)+(alfa*Pr*E*x*dc)-(((DT*Bs*Fiinf)/(Tinf*Dfi))*(dt+dy(4))))-(dc);
dy(7) = -1*(4*E+df+f);
dy(8) = (4*E*E)+((E*x+f-1)*(-1*(4*E+df+f)));
dy(9)=dg; %is it coorect for the definition of g'??
darova
darova 2020 年 4 月 14 日
But what is g? Looks like it's constant. If it exists the previous version was correct also
T K
T K 2020 年 4 月 14 日
The variable g is not constant. But the question about definition g 'in its equation f'' and y (9) is
T K
T K 2020 年 4 月 14 日
And dy(9) Please help me Dr darova??
darova
darova 2020 年 4 月 14 日
It's ok i think. I tried ode45 to see how results can look like
[x,y] = ode45(@projfun,[0 1],y0);
plot(x,y)
str = strsplit(sprintf('y%d ',1:9),' ');
legend(str)
I did my best. Don't know how proceed
T K
T K 2020 年 4 月 14 日
Sorry for your effort, Doctor darova Finally, you can send me the program with the ODE45 that I mentioned in the drawing
darova
darova 2020 年 4 月 14 日
I tried to reduce timespan
m = linspace(0,1);
T K
T K 2020 年 4 月 14 日
Finally, please doctor darova The code has been successfully created and given ode45 format and bvp4c function ??? First solution y0 = [1 0 1 1 1 1 1 1 0]; % Pr = 4; % [x,y] = ode45(@projfun,[0 .1],y0); % plot(x,y) % str = strsplit(sprintf('y%d ',1:9),' '); % legend(str)
Second solution options = bvpset('stats','on','RelTol',1e-4); m = linspace(0,1); solinit = bvpinit(m,y0);
darova
darova 2020 年 4 月 14 日
darova
darova 2020 年 4 月 14 日
T K
T K 2020 年 4 月 14 日
Thank you very much Doctor Best vote
darova
darova 2020 年 4 月 14 日
Don't forget about thumb up
darova
darova 2020 年 4 月 18 日
How can i trust you?
T K
T K 2020 年 4 月 18 日
There is no problem and I am very happy for you

サインインしてコメントする。

カテゴリ

ヘルプ センター および File ExchangeMathematics についてさらに検索

質問済み:

T K
2020 年 4 月 13 日

コメント済み:

T K
2020 年 4 月 18 日

Community Treasure Hunt

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

Start Hunting!

Translated by