2
$\begingroup$

I've tried to create the following example. Suppose that I have the differential equation:

$$ U''(x) = \frac{U'(x) - U(x)}{x}, $$

which I know 1 boundary condition and know that this function should behave normal at the origin. Thus, I would assume that to keep the function normal it should be true that:

$$ U'(0) = U(0), $$

with say my 1 boundary condition is:

$$ U(1) = 10, $$

or something of the sort. I would have then expected I could solve this differential equation with:

soln = NDSolve[{U''[x] == (U'[x] - U[x])/x, U'[0] == U[0], 
   U[10] == 10}, 
  U, 
  {x, 0, 1},
  Method -> {"Shooting", 
         "StartingInitialConditions" -> {U[0] == 1}}]

however, this immediately gives an infinite expression. Looking at the trace it seems that it does not consider that the numerator is $+1-1$ at all:

(-1. + 1.)/0.
Message[Power::infy, 1/0.]

Any thoughts on how to analyze such a differential equation? Any response is greatly appreciated.

$\endgroup$

3 Answers 3

3
$\begingroup$

Since the deq is singular at zero, the easiest workaround is to start the solution a little above zero.

Assign an small value

\[Epsilon] = 10^-10

soln = NDSolve[{U''[x] == (U'[x] - U[x])/x, U[\[Epsilon]] == 1, 
    U[1] == 10}, U[x], {x, \[Epsilon], 1}, 
   Method -> {"Shooting", 
     "StartingInitialConditions" -> {U'[\[Epsilon]] == 1}}, 
   WorkingPrecision -> 25] // Flatten

I raised the working precision a little to get rid of some warnings about precision.

U[x_] = U[x] /. soln

This particular problem can also be solved analytically:

sol = DSolve[{deq, u'[0] == u[0], u[0] == 1, u[1] == 10}, 
   u[x], {x, 0, 1}] // Flatten

u[x_] = u[x] /. sol

Check

Limit[u'[x] - u[x], x -> 0]
(*  0  *)

u[1]
(*  10  *)

Limit[u[x], x -> 0]
(*  1  *)

Mathematica gets infinity for x = 0, but taking the limit yields correct values.

Plot the numeric and symbolic solutions.

Plot[{U[x], u[x]}, {x, 0, 1}]

enter image description here

And they overlay.

$\endgroup$
3
$\begingroup$

The ODE may be solved symbolically with

uFN = DSolveValue[{U''[x] == (U'[x] - U[x])/x}, U, {x, 0, 2}]

(* Function[{x}, 2 x (BesselJ[2, 2 Sqrt[x]] C[1] - BesselY[2, 2 Sqrt[x]] C[2])]  *)

The singularity enforces the desired "boundary condition" at x == 0 on the whole solution space:

Limit[uFN'[x] - uFN[x], x -> 0, Direction -> "FromAbove"]
(*  0  *)

The BC at x == 10 reduces the dimension of the solution space by 1:

uFN10 = DSolveValue[{U''[x] == (U'[x] - U[x])/x, U[1] == 10}, 
  U, {x, 0, 2}]
(*
Function[{x}, 
 2 x (BesselJ[2, 2 Sqrt[x]] C[1] + 
   (BesselY[2, 2 Sqrt[x]] (5 - BesselJ[2, 2] C[1]))/BesselY[2, 2])]
*)

Plot[Table[uFN10[x] /. C[1] -> c, {c, 0, 10, 2}] // Evaluate, {x, 0, 1.1}]

enter image description here

$\endgroup$
1
  • 1
    $\begingroup$ BTW, to see the family of solutions with NDSolve, just use an IVP at U[10] == 10 (or u[1] == 10 as above), with a varying value for U'[10]. This make an nice, suggestive image: ListLinePlot@Table[NDSolveValue[{U''[x] == (U'[x] - U[x])/x, U'[10] == m, U[10] == 10}, U, {x, 0 + $MachineEpsilon, 25}], {m, -3, 3}] $\endgroup$
    – Michael E2
    Dec 7, 2021 at 2:59
1
$\begingroup$

Start NDSolve at a slightly negative x with m as parameter for U'[1] == m , and you have no problems at x == 0. Determine the corresponding m at x==0 for U[m][0]== "anyValue".

As @MichaelE2 says, "The singularity enforces the desired "boundary condition" at x == 0"

sol = Solve[U''[x] == (U'[x] - U[x])/x, U[x]]
sol /. x -> 0
(*   {{U[0] -> Derivative[1][U][0]}}   *)

Usol[m_?NumericQ] := 
U /. First@
NDSolve[{U''[x] == (U'[x] - U[x])/x, U'[1] == m, U[1] == 10}, 
U, {x, -10^-14, 3/2}]

Plot[Evaluate[Table[Usol[m][x], {m, -50, 50, 11/3}]], {x, 0, 1.5}]

enter image description here

Max@Abs@Table[Usol[m][0] - Usol[m]'[0], {m, -50, 50, 11/3}]
(*   2.30926*10^-13   *)

m1 = m /. FindRoot[Evaluate[Usol[m][0] == 1], {m, 3}]
(*   13.5113   *)

m1000 = m /. FindRoot[Evaluate[Usol[m][0] == 1000], {m, 3}]
(*   -2817.85   *)

Plot[Evaluate[Usol[m1000][x]], {x, 0, 1}]
$\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.