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$
    – user49048
    Commented Mar 11, 2022 at 17:47
  • 1
    $\begingroup$ I am using version 13.0 $\endgroup$
    – Johnson
    Commented Mar 11, 2022 at 17:59
  • $\begingroup$ thanks for that! $\endgroup$
    – user49048
    Commented Mar 11, 2022 at 17:59
  • $\begingroup$ @kcr Does the solution code from this work? $\endgroup$
    – Johnson
    Commented Mar 11, 2022 at 18:03
  • $\begingroup$ yup, the answer from that link works fine on v12 $\endgroup$
    – user49048
    Commented Mar 11, 2022 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$
    – Johnson
    Commented Mar 15, 2022 at 15:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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