how to solve this bvp4c?

Question

NURUL AISYAH JOHAN 2020-1-5

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/498953-how-to-solve-this-bvp4c

评论： Didarul Ahasan Redwan 2020-4-7

采纳的回答： darova

E.m

在 MATLAB Online 中打开

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