2
$\begingroup$

My code is working fine until I try to increase the value of parameter phi1, error such as

NDSolveValue::ndsz: At x == 0.23372152690117168`, step size is effectively zero; singularity or stiff system suspected.

starts to appear. Previously, I fixed this error in my code by employing this trick. However, this does not work anymore when i increase the value of parameter phi1

Clear["Global`*"]
(*constants*)
phi1 = 12;
KcKbRatio = 2;
DbDaRatio = 2;
DcDaRatio = 0.5;
Cref = 1;
lambda = 2;
aEnd = 0.95;
bEnd = 0.8;
cEnd = 0.6;
alpha = 1 - cEnd *Cref/lambda - 
   KcKbRatio*Cref/(KcKbRatio + 1)/DcDaRatio/lambda*aEnd ;
del = $MachineEpsilon; 
(*sidenotes:
when lambda>Cref, graphs tend to look normal
*)

(*the system of ode and bcs*)
ode = {a''[x] + 2/x*a'[x] - phi1^2*(1 + KcKbRatio)*phi[x]*a[x] == 0,
    b''[x] + 2/x*b'[x] + phi1^2/DbDaRatio*phi[x]*a[x] == 0,
   phi''[x] + 2/x*phi'[x] - 
     KcKbRatio*phi1^2*Cref/DcDaRatio/lambda*phi[x]*a[x] == 0};

bcs = {a'[del] == 0, a[1] == aEnd , b'[del] == 0, b[1] == bEnd, 
   phi'[del] == 0, phi[1] == (lambda - cEnd *Cref)/lambda};

(*ndsolve*)
{asol, bsol, phisol} = 
  NDSolveValue[{ode, bcs}, {a, b, phi}, {x, del, 1}];

(*Dsolve*)
exactsol =
  DSolve[{a''[x] + 2/x*a'[x] - 
      phi1^2*(1 + KcKbRatio)*(1 - cEnd*Cref/lambda)*a[x] == 0, 
    a'[0] == 0, a[1] == aEnd, 
    b''[x] + 2/x*b'[x] + 
      phi1^2/DbDaRatio*(1 - cEnd*Cref/lambda)*a[x] == 0, b'[0] == 0, 
    b[1] == bEnd}, {a[x], b[x]}, x];

(*Plot*)
Plot[{asol[x], bsol[x], phisol[x], 
  KcKbRatio*Cref/(KcKbRatio + 1)/DcDaRatio/lambda*asol[x] + alpha, 
  Evaluate[{a[x], b[x]} /. exactsol]}, {x, del, 1},
 PlotLegends -> "Expressions"]
$\endgroup$
7
  • $\begingroup$ I quickly tried the code you provided on version 12, and I get the stifness error even for phi1 = 0.01. Is it possible that you missed some part of your code? $\endgroup$
    – kcr
    Mar 11 at 17:47
  • 1
    $\begingroup$ I am using version 13.0 $\endgroup$
    – joe
    Mar 11 at 17:59
  • $\begingroup$ thanks for that! $\endgroup$
    – kcr
    Mar 11 at 17:59
  • $\begingroup$ @kcr Does the solution code from this work? $\endgroup$
    – joe
    Mar 11 at 18:03
  • $\begingroup$ yup, the answer from that link works fine on v12 $\endgroup$
    – kcr
    Mar 11 at 18:04

1 Answer 1

2
$\begingroup$

It works when using a piecwise expression for x -> 0.

Since a'[0] == 0 the expression 2 a'[x]/x tends to 1. Therefore substitute for a and equivalent for b and phi

(2 Derivative[1][a][x])/x -> 
     Piecewise[{{1, x < $MachineEpsilon}}, (2 Derivative[1][a][x])/x];

Clear["Global`*"]
(*constants*)
phi1 = 12;
KcKbRatio = 2;
DbDaRatio = 2;
DcDaRatio = 0.5;
Cref = 1;
lambda = 2;
aEnd = 0.95;
bEnd = 0.8;
cEnd = 0.6;
alpha = 1 - cEnd*Cref/lambda - 
   KcKbRatio*Cref/(KcKbRatio + 1)/DcDaRatio/lambda*aEnd;
del = $MachineEpsilon;
(*sidenotes:when lambda>Cref,graphs tend to look normal*)

(*the system of ode and bcs*)
ode = {a''[x] + 2/x*a'[x] - phi1^2*(1 + KcKbRatio)*phi[x]*a[x] == 0, 
   b''[x] + 2/x*b'[x] + phi1^2/DbDaRatio*phi[x]*a[x] == 0, 
   phi''[x] + 2/x*phi'[x] - 
     KcKbRatio*phi1^2*Cref/DcDaRatio/lambda*phi[x]*a[x] == 0};

bcs = {a'[del] == 0, a[1] == aEnd, b'[del] == 0, b[1] == bEnd, 
   phi'[del] == 0, phi[1] == (lambda - cEnd*Cref)/lambda};

   ode2 = {-432*a[x]*phi[x] + Piecewise[
       {{1, x < $MachineEpsilon}}, (2*Derivative[1][a][x])/x] + 
      Derivative[2][a][x] == 0, 
    72*a[x]*phi[x] + Piecewise[{{1, x < $MachineEpsilon}}, 
       (2*Derivative[1][b][x])/x] + Derivative[2][b][x] == 0, 
    -288*a[x]*phi[x] + Piecewise[{{1, x < $MachineEpsilon}}, 
       (2*Derivative[1][phi][x])/x] + Derivative[2][phi][x] == 
     0}; 

Method "StiffnessSwitching" produced quite large errors for lhs of ode. Method "BDF" performes better, as seen by the second graph.

ndsol = NDSolve[Join[ode2, bcs], {a, b, phi}, {x, 0, 1}, 
          Method -> "BDF"]

Plot[Evaluate[{a[x], b[x], phi[x]} /. ndsol[[1]]], {x, 0, 1}, 
 PlotStyle -> {Red, Green, Blue}, PlotRange -> All, 
 GridLines -> {Automatic, {19/20, 4/5, 7/10}}]

enter image description here

Plot[Evaluate[ode[[All, 1]] /. ndsol[[1]]], {x, 0, 1}, 
  PlotStyle -> {Red, Green, Blue}, PlotRange -> .02]

enter image description here

May be a more smooth approach to zero than piecwise could improve solutions near zero.

.

$\endgroup$
1
  • $\begingroup$ It is not working for my mathematica 13.0.1 for Microsoft Windows (64-bit) The following error would occur when i copy this into my mathematica: $\endgroup$
    – joe
    Mar 15 at 15:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.