Is my code for bvp4c for a set of odes is correctly written
3 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
I have the following set of equations:
f"=g(g^2+gamma^2)/(g^2+lambda*gamma^2) ------ (1)
g'= (1/3)*f'^2-(2/3)*(f*f")+ Mn*f' ------------------------(2)
t"+Rd*t"+ 2*Pr*f*t'/3+ Nb*t'*p'+Nt*(t')^2= 0------(3)
p"+(2*Lew*f*p')/3+ Nt*t"/Nb= 0 ------------------------(4)
'f', 'g', 't', 'p' are the dependent variables and η is the independent variable
And the boundary conditions are as follows:
f=0, f'= delta, t=1, p=1 at y=0
f'-->0, t-->0, p-->0 as y----> infinity
So solving these equations numerically I have used 'bvp4c' solver
So converting these higher order equations to first order equations, I have let
f = y1, df/dη = y2, g = y3, t = y4, dt/dη = y5, p = y6, dp/dη = y7
So the first order equations are
y1' = y2,
y2' = {y3*(y3^2+ gamma^2)}/(y3^2+ lambda*gamma^2) ,
y3' = y2^2/3 − (2/3)*y1*{y3*(y3^2+ gamma^2)/(y3^2+ lambda*gamma^2)} + Mny2,
y4' = y5,
y5' = −(1/(1+Rd))*{(2*Pr*y1*y5)/3 + Nb*y5*y7 + Nt*y5^2} ,
y6' = y7,
y7' = −(2*Le*y1*y7)/3 + {Nt/Nb*(1+Rd)}{(2*Pr*y1*y5)/3 + Nb*y5*y7 + Nt*y5^2}
I would like to know if my written code is correct and if I have got the correct plots
function sol= proj
clc;clf;clear;
global lambda gama Pr Rd Lew Nb Nt Mn m
gama=1;
Lew=1;
Nt=1;
Rd=1;
Pr=5;
Mn=1;
pp=[0.5 1 1.5];
qq=[0.2 0.5 0.9];
%figure(1)
%plot(2,1);hold on
options=bvpset('stats','on','RelTol',1e-9);
m=linspace(0,10);
solinit= bvpinit(m,[1,0,0,0,0,0,0]);
for i=1:numel(pp)
lambda=pp(i);
for i=1:numel(qq);
Nb=qq(i)
sol= bvp4c(@projfun,@projbc,solinit,options);
y1=deval(sol,0)
solinit= sol;
plot(sol.x,sol.y(6,:));hold on
end
end
end
function f= projfun(x,y)
global lambda gama Pr Rd Lew Nb Nt Mn
f= [y(2)
y(3)*(y(3)^2+gama^2)/(y(3)^2+lambda*gama^2)
y(2)^2/3-(2*y(1)*y(3)*(y(3)^2+gama^2))/(3*(y(3)^2+lambda*gama^2))+Mn*y(2)
y(5)
-(2*Pr*y(1)*y(5))/(3*(1+Rd)) - (Nb*y(5)*y(7))/(1+Rd) - (Nt*y(5)^2)/(1+Rd)
y(7)
-(2*Lew*y(1)*y(7))/3+ Nt*((2*Pr*y(1)*y(5))/(3*(1+Rd)) + (Nb*y(5)*y(7))/(1+Rd) + (Nt*y(5)^2)/(1+Rd))/Nb];
end
function res= projbc(ya,yb)
res= [ya(1); ya(2)-1; ya(4)-1.0; ya(6)-1.0; yb(2); yb(4); yb(6)];
end
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/171892/image.jpeg)
2 comentarios
Walter Roberson
el 9 de Abr. de 2018
If you have the Symbolic Toolbox, then I recommend that you look at odeFunction and go through the series of steps shown in the documentation to convert an ode right through to a function handle that can be used with the ode* routines.
Mayokun Ojediran
el 3 de Jul. de 2019
Hello, I am trying to adapt your code to solve the following equations:
(2n+1/2(n+1))f*theta' = (1/Bo^(2/n+1))*theta''
1/2*f'^2 +(2n+1/2(n+1))*f*f'' = -(f''^(n-1)*f'' + f''^(n-1)*f''') + theta
'f' and 'theta' are the dependent variables and η is the independent variable and Bo is a constant
f= f' = 0 , theta=1 at y=0
f'-->0, theta-->0 as y----> infinity
And the boundary conditions are as follows:
Converting to First order equations,
theta = y1, dtheta/dη = y2, f = y3, df/dη = y4, d^2f/d\η^2 = y5
The first order equations are
y(2);
Bo^(2/(n+1))*((2*n+1)/2*(n+1))*((2*n+1)/2*(n+1))*y(3)*y(2);
y(4);
y(5);
(-y(4)^2*n+2*y(3)*y(5)*n+2*y(1)*n-y(4)^2+y(3)*y(5)+2*y(1))/(2*y(5)^(n-1)*n*(n+1));
I have tried adapting your code but not completely sure how to implement the boundary conditions, please help.
function sol= proj
clc;clf;clear;
global n Bo
n = 0.4;
Bo = 0.0003548133892;
pp=[0.5 1 1.5];
qq=[0.2 0.5 0.9];
%figure(1)
%plot(2,1);hold on
options=bvpset('stats','on','RelTol',1e-9);
m=linspace(0,10);
solinit= bvpinit(m,[1,0,0,0,0,0,0]);
for i=1:numel(pp)
lambda=pp(i);
for i=1:numel(qq);
Nb=qq(i)
sol= bvp4c(@projfun,@projbc,solinit,options);
y1=deval(sol,0)
solinit= sol;
plot(sol.x,sol.y(4,:));hold on
end
end
end
function f= projfun(x,y)
global lambda gama Pr Rd Lew Nb Nt Mn
f= [y(2)
Bo^(2/(n+1))*((2*n+1)/2*(n+1))*((2*n+1)/2*(n+1))*y(3)*y(2);
y(4)
y(5)
(-y(4)^2*n+2*y(3)*y(5)*n+2*y(1)*n-y(4)^2+y(3)*y(5)+2*y(1))/(2*y(5)^(n-1)*n*(n+1))];
end
function res= projbc(ya,yb)
res= [ya(1); ya(2)-1; ya(4)-1.0; ya(6)-1.0; yb(2); yb(4); yb(6)];
end
Respuestas (0)
Ver también
Categorías
Más información sobre Ordinary Differential Equations en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!