Group Abstract

Message Boards

WOLFRAM COMMUNITY

5.2K Views

3 Replies

1 Total Like

View groups...

Follow this post

Share this post:

GROUPS:

Solve a hypergeometric ODE NDSolve?

Brad Klee

Posted 6 years ago

POSTED BY: Brad Klee

3 Replies

Sort By:

Brad Klee

Posted 6 years ago

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]] 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] Compare with earlier error plot. This quick hack is already doing better than NDSolve, and accepts initial condition at x=0. Any thoughts?

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

Brad Klee

Posted 6 years ago

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]]}] 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}] 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

Gianluca Gorni

Gianluca Gorni, University of Udine

Posted 6 years ago

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

Reply Preview

Attachments

Remove Add a file to this post

Follow this discussion

or Discard

Feedback