How to make NDSolve solve this apparently easy DE system?

Questions
, ,
1

Hello,

I have been trying to solve a BVP for a system of two apparently easy DEs, by using the code:

k=7;   (OK)
(k=8; fails )g=1;
d=1;
sol=NDSolve[{u''[x]-ku[x]/(g+u[x])==0,v''[x]+ku[x]/(g+u[x])/d==0,u'[0]==0,u[1]==1,v[0]==0,v[1]==0},{u,v},{x,0,1},
WorkingPrecision->80,PrecisionGoal->20,MaxSteps->∞];
ufun[x_]:=Evaluate[u[x]/.sol];
vfun[x_]:=Evaluate[v[x]/.sol];
Plot[{ufun[x],vfun[x]},{x,0,1}]

I need to solve the system for large k values, reaching, say 10^4 or higher, with a precision of at least 20 digits. Unfortunately, this code works for 0<=k<=7, but fails to produce results for k>=8. Looking at the solutions obtained for k=7, they don’t seem particularly challenging, except for the fact that v[0]=v[1]=0, so that the precision requested cannot be obviously achieved at x=0 or 1. But error message suggests some convergence problem in FindRoot. My question is: what would be the best choice of the method, and parameters, to achieve my goals, and why the current choice does not work?

Lesław

2

I have no problems with this:

k = 8; g = 1; d = 1;
{ufun, vfun} = 
  NDSolveValue[{u''[x] - k u[x]/(g + u[x]) == 0, 
    v''[x] + k u[x]/(g + u[x])/d == 0,
    u'[0] == 0, u[1] == 1, v[0] == 0, v[1] == 0},
   {u, v}, {x, 0, 1}];
Plot[{ufun[x], vfun[x]}, {x, 0, 1}]
3

Yes, but how to make it effective to provide 20 accurate digits that I need?

Leslaw

4
Open in Wolfram Cloud ↗
5

You can concentrate on u[x], becausev[x] can be deduced from u[x] with

v[x] == (u[0] - x u[0] + x u[1] - u[x])/d

Unfortunately, I cannot help you with u[x]

k = 8; g = 1;
ufun = NDSolveValue[{u''[x] - k u[x]/(g + u[x]) == 0,
   u'[0] == 0, u[1] == 1}, u, {x, 0, 1},
  WorkingPrecision -> 300, PrecisionGoal -> 20, MaxSteps -> Infinity]
{ufun'[0], ufun[1]}
6

You have commented out one boundary condition: u’[0]==0. How can this lead to a correct solution?

Lesław

7

Gianluca: How did you obtain this relationship:

Obviously , u[x]+d*v[x] must be a linear function of x, say u[x]+d*v[x] ==A x + B, and then from the BCs at x=1 we get B==1-A. But I don’t see a way to obtain A.

Lesław

8

One way is this

DSolve[{u''[x] + d v''[x] == 0, v[0] == 0, v[1] == 0},
 v[x], x]

otherwise

u[x] + d*v[x] == A x + B /.
  {{x -> 0}, {x -> 1}} /.
 {v[0] -> 0, v[1] -> 0}
Solve[%, {A, B}]
9

I have made additional tests, focusing on solving u[x], and considering a linearised problem obtained by omitting u[x] in the denominator of k*u[x]/(g + u[x]).

In such a case an analytical solution is obtainable, which allows one to check the actual relative error of the solution. The following code easily produces relative errors of u[x] smaller than 10^-20, even for g=10^4, as the last plot reveals. These errors depend on the choice of working precision, difference order, and step size. Additionally specifying precision goal and accuracy goal does not seem to have an effect.

uanal[x_]:=(E^(-Sqrt[(k/g)] (-1+x)) (1+E^(2 Sqrt[k/g] x)))/(1+E^(2 Sqrt[k/g]));

wp=200;
do=20;
ss=1/100;
k=10000;

g=1;
d=1;
sol=NDSolve[{u''[x]-k*u[x]/g==0,u'[0]==0,u[1]==1},u[x],{x,0,1},
WorkingPrecision->wp,
Method->{"FixedStep",Method->{"ImplicitRungeKutta","DifferenceOrder"->do}},StartingStepSize->ss,InterpolationOrder->All];
ufun[x_]:=Evaluate[u[x]/.sol];
vfun[x_]:=(ufun[0]+(ufun[1]-ufun[0])*x-ufun[x])/d;
Plot[{ufun[x],vfun[x]},{x,0,1},PlotRange->{0,1},PlotLegends->"Expressions"]
Plot[Log[10,Abs[ufun[x]/uanal[x]-1]],{x,0,1},WorkingPrecision->wp,PlotLegends->"Expressions"]

Unfortunately, the complete nonlinear DE makes the above code to freeze. Or, by modifying the parameter settings I get various error messages. Isn’t there any way to overcome these difficulties?

What is actually the algorithm used? Some sort of a shooting method?

Lesław