I want a code for solving a coupled 3rd order and 2nd order ode using shooting method and RK-4 numerical technique , please if anyone could help

3 ビュー (過去 30 日間)
(1+2M*eta)f''' + 2M*f"+ f*f" -f'^2 - k1*f' + lambda*theta=0 ---------- (1)
(1+2M*eta)theta" + 2M*theta' + Pr(f*theta'-f'*theta)=0 -------(2)
'f' and 'theta' are functions of 'eta', eta is an independent variable
3 initial conditions are given: eta=0, f(0)=0, f'(0)=1,theta(0)=1
Say I reduce these equations (1) and (2) to five ode (shooting method)
f'=z ; f(0)=0 -----(3)
z'=p ; z(0)=1 ------(4)
p'= (-2M*f"-f*f"+f'^2+k1*f'-lamda*theta)/(1+2*M*eta) ;p(0)= (guess value) -----(5)
theta'= q ; theta(0)=1 -----(6)
q' = (-2M*theta'-Pr(f*theta'-f'*theta))/(1+2*M*eta) ; q(0)= (guess value) ------(7)
The boundary conditions that needs to be satisfied are: f'(eta=10)= 0 and theta(eta=10)=0 as eta=10
Given:
M= 1
k1= 0.1
lamda= 0.1
Pr= 0.7
taking step length: h= 0.01
  4 件のコメント
Torsten
Torsten 2017 年 11 月 16 日
Here is the link to a very similar problem set up for BVP4C:
https://de.mathworks.com/matlabcentral/answers/366666-bvp4c-error-cannot-solve-collocation-equations-singular-jacobian
Best wishes
Torsten.
naygarp
naygarp 2017 年 11 月 21 日
Hello,
I have copied the code given on this link
https://de.mathworks.com/matlabcentral/answers/366666-bvp4c-error-cannot-solve-collocation-equations-singular-jacobian
I could not fully grasp which thing to modify apart from the function handle. Please could you show me where to put the initial conditions and how to get the results, I need two guess initial values for y(3) and y(5)..I have encountered a large no. of errors
global Pr k1 M lambda
Pr=0.7; k1=0.1; M=1; lambda=0.1;
%
rlow=0;
rhigh=10;
N=1000;
options=bvpset('stats','on','RelTol',1e-5);
solinit=bvpinit(linspace(rlow,rhigh,N),[0, -1, 0, 1 ,0]); %
sol=bvp4c(@projode,@mybcs,solinit,options);
function dy=projode(n,y)
global Pr k1 M lambda
dy=[y(2); y(3); (-2*M*y(3)-y(1)*y(3)+y(2)^2+k1*y(2)-lambda*y(4))/(1+2*M*n); y(5);... (-2*M*y(5)-Pr*f(1)*y(5)+Pr*y(2)*y(4))/(1+2*M*n)];
end
function res=mybcs(ya,yb)
res=[ya(1) ya(2) ya(4) yb(2) yb(4)];
end

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

採用された回答

Torsten
Torsten 2017 年 11 月 22 日
Try
function main
global Pr k1 M lambda
Pr=0.7; k1=0.1; M=1; lambda=0.1;
rlow=0;
rhigh=10;
N=1000;
options=bvpset('stats','on','RelTol',1e-5);
solinit=bvpinit(linspace(rlow,rhigh,N),[0, -1, 0, 1 ,0]); %
sol=bvp4c(@projode,@mybcs,solinit,options);
function dy = projode(n,y)
global Pr k1 M lambda
dy = [y(2); y(3); (-2*M*y(3)-y(1)*y(3)+y(2)^2+k1*y(2)-lambda*y(4))/(1+2*M*n); y(5); (-2*M*y(5)-Pr*y(1)*y(5)+Pr*y(2)*y(4))/(1+2*M*n)];
function res = mybcs(ya,yb)
res = [ya(1); ya(2)-1.0; ya(4)-1.0 ; yb(2); yb(4)];
Best wishes
Torsten.
  6 件のコメント
naygarp
naygarp 2017 年 11 月 28 日
Thanks, so I think the code is right if you confirm it...
Now I have run the code for various values of M from 0 to 1 (0,0.25,0.5,0.75,1)
And plotted the graphs of 'f' vs eta' and 'theta vs eta', the graphs I got and the graphs published in research papers are different.
Could you figure it out why
|function sol= proj
global M k1 p Pr
k1=0.1;
p=0.1;
Pr=0.7;
MM=[0:0.25:1];
figure(1)
subplot(2,1,1);
subplot(2,1,2);
solinit= bvpinit([0:0.01:10],[0,1,0,1,0]);
for M= MM
sol= bvp4c(@projfun,@projbc,solinit);
solinit= sol;
subplot(2,1,1);plot(sol.x,sol.y(2,:));hold on
subplot(2,1,2);plot(sol.x,sol.y(4,:));hold on
end
end
function f= projfun(x,y)
global M k1 p Pr
f= [y(2);y(3);(y(2)^2+k1*y(2)-2*M*y(3)-p*y(4))/(1+2*M*x);y(5);(Pr*y(2)*y(4)-Pr*y(1)*y(5)-2*M*y(5))/(1+2*M*x)];
end
function res= projbc(ya,yb)
res= [ya(1); ya(2)-1.0; ya(4)-1.0 ; yb(2); yb(4)];
end|
Torsten
Torsten 2017 年 11 月 28 日
If you include the research paper and your graphs: maybe.
Best wishes
Torsten.

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

その他の回答 (2 件)

naygarp
naygarp 2017 年 11 月 28 日
  7 件のコメント
naygarp
naygarp 2017 年 11 月 30 日
Well thanks again for the confirmation, I would like to ask again
Do you think adding this line in the code
'options=bvpset('stats','on','RelTol',1e-9)'
and putting it in
'sol=bvp4c(@projode,@mybcs,solinit,options)'
would do any corrections.
Torsten
Torsten 2017 年 12 月 1 日
I think no, but don't you have MATLAB available to test it ?
Best wishes
Torsten.

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


iqra
iqra 2024 年 1 月 31 日
function main
global Pr k1 M lambda
Pr=0.7; k1=0.1; M=1; lambda=0.1;
rlow=0;
rhigh=10;
N=1000;
options=bvpset('stats','on','RelTol',1e-5);
solinit=bvpinit(linspace(rlow,rhigh,N),[0, -1, 0, 1 ,0]); %
sol=bvp4c(@projode,@mybcs,solinit,options);
function dy = projode(n,y)
global Pr k1 M lambda
dy = [y(2); y(3); (-2*M*y(3)-y(1)*y(3)+y(2)^2+k1*y(2)-lambda*y(4))/(1+2*M*n); y(5); (-2*M*y(5)-Pr*y(1)*y(5)+Pr*y(2)*y(4))/(1+2*M*n)];
function res = mybcs(ya,yb)
res = [ya(1); ya(2)-1.0; ya(4)-1.0 ; yb(2); yb(4)];

カテゴリ

Help Center および File ExchangeNumerical Integration and Differential Equations についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by