1
$\begingroup$

I am new to mathematica and I am trying to solve the systems of ODE as shown below

Clear["Global`*"]
(*constants*)
phi1 = 1;
s = 0.5;
P = 1;
Cref = 1;
Q = 1;
lambda = 1;

(*the system of ode*)
ode = {x*a''[x] - 2*a'[x] - phi1^2*(1 + s)*phi[x]*a[x] *x == 0,
    x*b''[x] - 2*b'[x] - phi1^2/P*phi[x]*a[x]*x == 0,
   x*phi''[x] - 2*phi'[x] - s*phi1^2*Cref/Q/lambda*phi[x]*a[x]*x == 0};
bcs = {a'[0] == 0, a[1] == 0.9, b'[0] == 0, b[1] == 0.95, 
   phi'[0] == 0, phi[1] == (lambda - 0.6*Cref)/lambda};

(*ndsolve*)
{asol, bsol, phisol} = 
  NDSolveValue[{ode, bcs}, {a, b, phi}, {x, 0, 1}];
Method -> {"MethodOfLines", 
   "SpatialDiscretization" -> {"FiniteElement"}};
Plot[{asol[x], bsol[x], phisol[x]}, {x, 0, 1}, 
 PlotLegends -> "Expressions"]

I would get error 'Power::infy' and 'Infinity::indet' when i am solving the system from x = 0 to 1. However, this can be overcome when i replace 0 with a close to 0 value. May i ask if there is other better solution for this?

Power::infy: Infinite expression 1/0. encountered.

Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered.

$\endgroup$

1 Answer 1

1
$\begingroup$

When this happens means you have issue with starting at zero, where first step gives 1/0 (since you are starting from x=0).

Without looking more at it to find why, and which term, a quick workaround is to start from little after x=0. Using $MachineEpsilon actually worked.

Clear["Global`*"]
(*constants*)
phi1 = 1;
s = 0.5;
P = 1;
Cref = 1;
Q = 1;
lambda = 1;
del = $MachineEpsilon; (*added this *)

(*the system of ode*)
ode1 = x*a''[x] - 2*a'[x] - phi1^2*(1 + s)*phi[x]*a[x]*x == 0
ode2 = x*b''[x] - 2*b'[x] - phi1^2/P*phi[x]*a[x]*x == 0
ode3 = x*phi''[x] - 2*phi'[x] - s*phi1^2*Cref/Q/lambda*phi[x]*a[x]*x ==
   0
bcs = {a'[del] == 0, a[1] == 0.9, b'[del] == 0, b[1] == 0.95, 
  phi'[del] == 0, phi[1] == (lambda - 0.6*Cref)/lambda}

(*ndsolve*)
NDSolveValue[{ode1, ode2, ode3, bcs}, {a, b, phi}, {x, del, 1}]

Mathematica graphics

No errors now.

Again, the issue is that there is a term or derivative in your equations which gives 1/x when x=0 this gives a problem. So more investigation is need to find it.

$\endgroup$

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.