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

(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

[Torsten]

You are not allowed to use a numerical method different from RK4 in combination with shooting ?

If you are free in the choice of the numerical method, use MATLAB's BVP4C (as already suggested). The length of your code for the above problem will be about 10 lines. Just make an attempt.

Best wishes

Torsten.

[naygarp]

I am told to use only RK4 method with shooting as it is mostly followed in research papers. Could you please show me some lines on how to use BVP4C on this problem.

[Torsten]

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]

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

Answer 1

[Torsten]

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.

[naygarp]

In the line:

solinit=bvpinit(linspace(rlow,rhigh,N),[0, -1, 0, 1 ,0]);

what does the matrix [0, -1, 0, 1,0] mean?

my initial conditions are n=0, y(1)=0, y(2)=1, y(3)= guess value, y(4)=0, y(5)=1, y(6)= guess value.

And I want the final boundary condition at n=10, to be y(2)=0 and y(5)= 0

How should I make the adjustments

[Torsten]

I only see five equations in your code - so what is y(6) for which you want to prescribe an initial guess ?

In general,

[0, -1, 0, 1, 0]

in

solinit=bvpinit(linspace(rlow,rhigh,N),[0, -1, 0, 1 ,0]);

means that you give the initial guesses for the solution variables for the complete interval [rlow rhigh] to be 0 for y(1), -1 for y(2), 0 for y(3), 1 for y(4) and 0 for y(5).

Additionally, if you want to set boundary values b1,...,b3 at r=0 for y(1),...,y(3) and boundary values b4,b5 at r=10 for y(4),y(5), e.g., then "res" in "mybcs" looks like

res=[ya(1)-b1,ya(2)-b2,ya(3)-b3,yb(4)-b4,yb(5)-b5]

Best wishes

Torsten.

[naygarp]

Here's a code that I managed to run without any errors

function sol= proj

global M k1 p Pr

M=0;

k1=0.1;

p=0.1;

Pr=0.7;

solinit= bvpinit([0:0.01:10],[0,1,0,1,0]);

sol= bvp4c(@projfun,@projbc,solinit)

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

Please check to see that I got the 'solinit' and 'res' vectors correct.

[Torsten]

This is the problem setup for

(1+2*M*eta)*f'''+2*M*f''+f*f''-f'^2-k1*f'+lambda*theta=0
(1+2*M*eta)*theta''+Pr*(f*theta'-f'*theta)+2*M*theta'=0

with boundary conditions

f(0)=0
f'(0)=1
theta(0)=1
f'(10)=0
theta(10)=0

and parameters

M=0
k1=0.1
lambda=0.1
Pr=0.7

If this is what you wanted to set, everything is fine.

Best wishes

Torsten.

[naygarp]

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]

If you include the research paper and your graphs: maybe.

Best wishes

Torsten.

Answer 2

[naygarp]

[naygarp]

[naygarp]

my graphs or you could run the code to get the graphs that I have provided with a modicification of including a 'for' loop-

[Torsten]

To check the equations, I need the research paper.

I solved your equations using a different numerical tool (COMSOL) and I get the same results as you do.

So my guess is that the authors of the paper solved equations that somehow differ from yours.

Best wishes

Torsten.

[naygarp]

I am attaching the research paper, please have a look.

[Torsten]

You forgot one term in the definition of f''', but this doesn't change the curves drastically:

f = [y(2);y(3);(y(2)^2+k1*y(2) -y(1)*y(3) -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)];

Since we got similar results from two different software programs, my guess is that the results in the article are wrong.

Best wishes

Torsten.

[naygarp]

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]

I think no, but don't you have MATLAB available to test it ?

Best wishes

Torsten.

Answer 3

[iqra]

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)];