How do I extract a particular value in a 'bvp4c' code solver?

3 visualizzazioni (ultimi 30 giorni)
I have this bvp4c code for the following 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)
With boundary conditions:
f=0, f'=1, t=1, p=1; at η=0
f'=0, t=0, p=0; at η-->infinity
I have used bvp4c solver to solve it numerically
function sol= proj
clc;clf;clear;
global lambda gama Pr Rd Lew Nb Nt Mn m
%Nb=0.5;
gama=1;
Mn=1;
Rd=0.1;
Pr=1;
Nb=0.3;
Lew=10;
%lambda=1.5;
pp=0.5;
qq=[0.05:0.01:0.5];
%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);
Nt=qq(i)
sol= bvp4c(@projfun,@projbc,solinit,options);
y1=deval(sol,0)
%y2=deval(y(6,:),0)
solinit= sol;
plot(sol.x,sol.y(6,:),'LineWidth',2);hold on
%legend('N_{b} = 0.1','N_{b} = 0.3','N_{b} = 0.5');
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/Nb)*((2*Pr*y(1)*y(5))/(3*(1+Rd)) + (Nb*y(5)*y(7))/(1+Rd) + (Nt*y(5)^2)/(1+Rd))];
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
I want to extract the values of y(6,:) particularly
Using 'deval(sol,0)' gives all the values of the 'y' rows
Is there any command to extract the values of y(6,:) or y(4,:) in particular.

Risposta accettata

MINATI
MINATI il 25 Apr 2018
Try with this f0 = deval(sol,0); f0(4,:)
Best Wishes
  1 Commento
naygarp
naygarp il 26 Apr 2018
Thanks, Pradyumna.
Could you please check, whether the above bvp4c code is correctly written, I mean the ODEfunction and the result obtained is correct.

Accedi per commentare.

Più risposte (0)

Prodotti

Community Treasure Hunt

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

Start Hunting!

Translated by