Message Boards Message Boards

Solve a hypergeometric ODE NDSolve?

Brad Klee, Wolfram Research
Posted 5 years ago
POSTED BY: Brad Klee
3 Replies
Sort By:

How To: Calculate All Derivatives

Mat = Outer[Coefficient,
    D[Dot[{1, D[-4 x (1 - x), x], -4 x (1 - x)},
        {T[x], D[T[x], x], D[T[x], {x, 2}]}],
       {x, #}] & /@ Range[0, 10], D[T[x],
       {x, #}] & /@ Range[0, 11]] /. x -> 0;
Prepend[Mat, Prepend[Table[0, 11], 1]] // MatrixForm
LinearSolve[Prepend[Mat, Prepend[Table[0, 11], 1]], 
 Prepend[Table[0, 11], 1]]

out

And notice, this acutally requires only one initial condition, not even two!

How To: DIY NDSolve

DerMat[n_] := Join[
   IdentityMatrix[n + 3][[1 ;; 2]],
   Outer[Coefficient,
    D[Dot[{1, D[-4 x (1 - x), x], -4 x (1 - x)},
        {T[x], D[T[x], x], D[T[x], {x, 2}]}],
       {x, #}] & /@ Range[0, n], D[T[x],
       {x, #}] & /@ Range[0, n + 2]]];

GetNext[mat_, cond_, dx_] := With[{Tx = Dot[
     LinearSolve[mat /. x -> cond[[3]], Join[cond[[1 ;; 2]],
       Table[0, Length[mat] - 2]]], 
     1/#! x^# & /@ Range[0, Length[mat] - 1]]
   }, {Tx, D[Tx, x], cond[[3]] + dx} /. x -> dx]

AbsoluteTiming[
 FApprox = Function[{MAT},
        Interpolation[
         NestList[GetNext[MAT, #, .0001] &, {1, 1/4, 0}, 
           9999][[All, {3, 1}]]]][DerMat[#]] & /@ {2, 3, 4, 5, 6} // 
    Quiet;]

Show[LogPlot[
    Abs[FApprox[[#]][x] - 2/Pi EllipticK[x]], {x, .001, .9999}, 
    PlotRange -> All, 
    PlotStyle -> {Red, Orange, Yellow, Green, Blue}[[#]]] & /@ 
  Range[5],
 LogPlot[10^(-#), {x, 0, 1}, PlotStyle -> Black] & /@ {7, 11},
 PlotRange -> All]

err2

Compare with earlier error plot. This quick hack is already doing better than NDSolve, and accepts initial condition at x=0. Any thoughts?
POSTED BY: Brad Klee

Your differential equation is of the second order, but you give three initial data. Also, at the the starting point x==0 of the equation is singular

The error message asked for "initial values for all the dependent variables ", then complained that I gave three I.C.s, even though they are consistent. You are right that x=0 is a singular point, but why should this matter at all? Look at the graphs:

Plot[{ LegendreQ[-(1/2), -1 + 2 x], EllipticK[x]}, {x, -1, 1}, 
 PlotStyle -> {Directive[Lighter@Lighter@Blue, Thick], 
   Directive[Thick, Black, Dashing[.05]]}]

EllipticK

This particular solution is smooth at x=0, so iterative solution algorithms should work fine. In fact, the function EllipticK has a known series expansion, so if nothing else, the NDSolve could expand the series and evaluate. Also compare with:

ICs = MapThread[Equal,
   {(D[T[x], {x, #}] /. x -> 10^(-10)) & /@  Range[0,1], 
    (D[(2/Pi) EllipticK[x], {x, #}] /. x -> 10^(-10)) & /@ 
     Range[0, 1]}];
F = T /. NDSolve[Prepend[ICs, Eq1], 
    T, {x, 10^(-10), 1 - 10^(-15)}][[1]];
LogPlot[F[x] - 2/Pi EllipticK[x], {x, 0, .99}]

error

Using approximately x=0 as the initial condition, we can compute an accurate but imprecise solution. So it looks like we are getting closer, but we still don't know:

  • Is it possible to set the initial condition at x=0 and calculate a numerical solution to arbitrary precision? If yes, how? If no, why not?

Unfortunately, this is starting to look more like the "Mission Critical" failure I was worried about earlier.
POSTED BY: Brad Klee

Your differential equation is of the second order, but you give three initial data. Also, at the the starting point x==0 of the equation is singular:

NDSolve[{T[x] + 4 (-1 + 2 x) Derivative[1][T][x] + 
    4 (-1 + x) x (T^\[Prime]\[Prime])[x] == 0, T[0] == 1, 
  T'[0] == 1/4, T''[0] == 9/32}, T, {x,
   0, 1}]

DSolve gives a solution, which is probably equivalent with yours:

DSolve[Prepend[ICs, Eq1], T, x]
POSTED BY: Gianluca Gorni
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract