how to solve this bvp4c?

조회 수: 3 (최근 30일)
NURUL AISYAH JOHAN
NURUL AISYAH JOHAN 2020년 1월 5일
댓글: Didarul Ahasan Redwan 2020년 4월 7일
Hello, can anyone help me solve this bvp4c program?
as to solve this problem, I already try numerous initial numbers.
it's either
Unable to solve the collocation equations -- a singular Jacobian encountered.
or
Warning: Unable to meet the tolerance without using more than 1250 mesh points.
The last mesh of 892 points and the solution are available in the output argument.
The maximum residual is 85.3834, while requested accuracy is 1e-07.
here is the bvp4c:
format long g
global alpha beta phi phid epsilon delta rhoS rhoF Y betaT Ec cpS cpF kf knf Pr Tv A B w Q ks
a=0; b=15;
phi=0.5;phid=0.2;delta=0.5;A=0.5;B=0.5;w=1.5;
betaT=0.2;alpha=0.5;beta=0.5;Ec=0.2;Y=0.01;Pr=6.2;Tv=0.488;epsilon=0.1;
rhoF=997.1;cpF=4179;kf=0.613;
rhoS=8933;cpS=385;knf=0.81629325;
Q=((ks+(2*kf))-(2*phi*(kf-ks)))/((ks+(2*kf))+(phi*(kf-ks))); ks=400;
solinit = bvpinit(linspace(a,b,100),@fluidparticle_init);
options = bvpset('stat','on','RelTol',1e-7);
sol = bvp4c(@fluidparticle_ode,@fluidparticle_bc,solinit,options);
sol.y(3,1)
sol.y(7,1)
plot(sol.x,sol.y(2,:),':r')
function dydx=fluidparticle_ode(x,y,alpha,beta,phi,phid,Y,betaT,Ec,kf,knf,Pr,Tv,A,B,rhoS,rhoF,cpS,cpF,Q)
global alpha beta phi phid Y betaT Ec kf knf Pr Tv A B rhoS rhoF cpS cpF Q
dydx=[y(2);
y(3);
(((1-phi)^2.5)*(1-phi+(phi*(rhoS/rhoF)))*((y(2)^2)-(y(1)*y(3))))-((((1-phi)^2.5)/(1-phid))*alpha*beta*(y(5)-y(2)));
y(5);
((y(5)^2)-(beta*(y(5)-y(2))))/y(4);
y(7);
(1-phi+(phi*((rhoS*cpS)/(rhoF*cpF))))*(Pr/Q)*((2*y(2)*y(6))-(y(1)*y(7)))-((Pr/Q)*((alpha*betaT*(y(8)-y(6)))-((alpha*Ec/Tv)*((y(5)-y(2))^2))))
((2*y(5)*y(8))+(Y*betaT*(y(8)-y(6))))/y(4)];
function res=fluidparticle_bc(ya,yb,epsilon,delta,w)
global epsilon delta w
res=[ya(2)-epsilon-(delta*ya(3))
ya(1)
yb(2)
yb(5)
yb(4)-yb(1)
ya(6)-1-(w*ya(7))
yb(6)
yb(8)];
function v=fluidparticle_init(x)
v=[0.5
-exp(-x)
exp(-x)
0.3
-0.5
-exp(-x)
exp(-x)
-exp(-x)];
what should I adjust in this program?
thank you in advance for your help and kindness.
  댓글 수: 2
darova
darova 2020년 1월 5일
Please attach the original equations
NURUL AISYAH JOHAN
NURUL AISYAH JOHAN 2020년 1월 6일
〖(1-ϕ_d)/(1-ϕ)^2.5〗f^''' - (1-ϕ_d) [1-ϕ+ϕ(ρ_s/ρ_f )](f^'2-〖ff〗^'' )+αβ(F^'- f^' )= 0.....(1)
F^'2-FF^''-β(f^'-F^' )=0.....(2)
(1/Pr)(k_nf/k_f ) θ^''-(1-ϕ+ϕ((ρC_p )_(s )/(ρC_p )_f ))(2f^' θ- fθ^' )+ αβ_T (θ_p-θ)+αEc/τ_v (F^'-f^' )^2 = 0.....(3)
2F^' θ_p-Fθ_p^'+γβ_T (θ_p-θ)=0.....(4)
Boundary condition:
f^' (0) = ε+δf^'' (0), f(0) = 0, θ(0) = 1+ωθ'(0)
f^' (η)=0, F^' (η)=0, F(η)=f(η), θ(η)=0, θ_p (η)=0 as η → ∞

댓글을 달려면 로그인하십시오.

채택된 답변

darova
darova 2020년 1월 6일
I re-wrote the main equations (compare yours and mine) and tried b = 1
Check the attachment
  댓글 수: 7
darova
darova 2020년 4월 7일
Why wasn't you? I don't know
I set up a limit b=6/7 and got this
Didarul Ahasan Redwan
Didarul Ahasan Redwan 2020년 4월 7일
can you please have a look at this?
I would be very thankful

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Logical에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by