       Re: NDSolve

• To: mathgroup at smc.vnet.net
• Subject: [mg21544] Re: NDSolve
• From: Paul Abbott <paul at physics.uwa.edu.au>
• Date: Fri, 14 Jan 2000 02:44:04 -0500 (EST)
• Organization: University of Western Australia
• References: <85esu5\$pca@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

```Klaus Wallmann wrote:

> Using MATHEMATICA Version 4.0 on a Pentium III PC with Windows NT Version
> 4, I tried to solve the following ordinary homogeneous non-linear second
> order differential equation with variable coefficients:
>
> y''[t] + p[t] y'[t] + q[t] y[t]/(y[t]+K)  = 0
>
> with boundary conditions
>
> y = 29
>
> y' = 0
>
> and
>
> K = 1
>
> p[t] = ((-1.15301 (-3 + ln[(0.769 + 0.102 exp(-0.036 t))^(2)]) + 0.2307
> exp(+0.036 t) (-1 + ln[(0.769 + 0.102 exp(-0.036 t))^(2)])^(2))/(314 (0.871
> + 0.769 (-1 + exp(+0.036 t) (-1 + ln[(0.769 + 0.102 exp(-0.036 t))^(2)]))
>
> q[t] = -[0.00203318 (0.231 - 0.102 exp(-0.036 t)) (0.102 (0.697676 -
> exp(-0.036 t))+ 0.008316 (93.2 + t))^(-1.142) (1 - ln[(0.769 + 0.102
> exp(-0.036 t))^(2)])]/(0.769 + 0.102 exp (-0.036 t))
>
> I tried NDSolve to solve this equation but I did not succeed.
>
> Is it possible to solve this equation or similar differential equations
> using the numerical procedures implemented in MATHEMATICA?

The following Notebook outlines one possible approach.

Cheers,
Paul

Notebook[{
Cell["\<\
Consider solving the following ordinary homogeneous non-linear \
second order differential equation with variable coefficients,\
\>", "Text"],

Cell[BoxData[
FormBox[
RowBox[{
RowBox[{"de", "=",
RowBox[{
RowBox[{\(\(\(y(t)\)\ \(q(t)\)\)\/\(y(t) + 1\)\), "+",
RowBox[{\(p(t)\), " ",
RowBox[{
SuperscriptBox["y", "\[Prime]",
MultilineFunction->None], "(", "t", ")"}]}], "+",
RowBox[{
SuperscriptBox["y", "\[Prime]\[Prime]",
MultilineFunction->None], "(", "t", ")"}]}], "==", "0"}]}],

Cell[TextData[{
"with boundary conditions ",
Cell[BoxData[
" and ",
Cell[BoxData[
FormBox[
RowBox[{
RowBox[{
SuperscriptBox["y", "\[Prime]",
MultilineFunction->None], "(", "6000", ")"}], "==", "0"}],
", where"
}], "Text"],

Cell[BoxData[
t_) = \(0.2307\ \[ExponentialE]\^\(\(+0.036\)\ t\)\ \((log(\((0.102\ \

\[ExponentialE]\^\(\(-0.036\)\ t\) + 0.769)\)\^2) - 1)\)\^2 - 1.15301\ \
\((log(\((0.102\ \[ExponentialE]\^\(\(-0.036\)\ t\) + 0.769)\)\^2) - 3)\)\)\/\

\(314\ \((0.769\ \((\[ExponentialE]\^\(\(+0.036\)\ t\)\ \((log(\((0.769 + \
0.102\ \[ExponentialE]\^\(\(-0.036\)\ t\))\)\^2) - 1)\) - 1)\) + 0.871)\)\); \

\)], "Input"],

Cell["and ", "Text"],

Cell[BoxData[
t_) = \(-\(\(0.00203318\ \((0.231 -
0.102\ \[ExponentialE]\^\(\(-0.036\)\ t\))\)\ \((1 -
log(\((0.769 + 0.102\ \[ExponentialE]\^\(\(-0.036\)\ t\))\)\

\^2))\)\)\/\(\((0.008316\ \((t + 93.2)\) + 0.102\ \((0.697676 - \
\[ExponentialE]\^\(\(-0.036\)\ t\))\))\)\^1.142\ \((0.769 +
0.102\ \[ExponentialE]\^\(\(-0.036\)\ t\))\)\)\)\); \)], \
"Input"],

Cell[TextData[{
Cell[BoxData[
FormBox[
StyleBox["NDSolve",
" cannot solve this boundary value problem directly. However it is possible
\
to solve such equations using ",
Cell[BoxData[
FormBox[
StyleBox["NDSolve",
"."
}], "Text"],

Cell[TextData[{
"Firstly, it is useful to graph the behaviour of ",
Cell[BoxData[
","
}], "Text"],

Cell[BoxData[
PlotRange -> All];\)\)], "Input"],

Cell[TextData[{
"and ",
Cell[BoxData[
":"
}], "Text"],

Cell[BoxData[
PlotRange -> All];\)\)], "Input"],

Cell[TextData[{
"We observe that ",
Cell[BoxData[
" is essentially flat for ",
Cell[BoxData[
"."
}], "Text"],

Cell[TextData[{
"Since ",
Cell[BoxData[
FormBox[
StyleBox["NDSolve",
" cannot solve this boundary value problem directly we use a shooting \
method instead. We can shoot forward (from the initial condition), "
}], "Text"],

Cell[BoxData[
FormBox[
RowBox[{\(forward[m_?NumberQ]\), ":=",
RowBox[{"y", "/.",
RowBox[{"First", "[",
RowBox[{"NDSolve", "[",
RowBox[{
RowBox[{"{",
RowBox[{"de", ",", \(y(10) == 29\), ",",
RowBox[{
RowBox[{
SuperscriptBox["y", "\[Prime]",
MultilineFunction->None], "(", "10", ")"}], "==",
"m"}]}], "}"}], ",", "y", ",", \({t, 10, 6000}\)}],

Cell[TextData[{
"supplying an initial slope, ",
Cell[BoxData[
", or backwards"
}], "Text"],

Cell[BoxData[
FormBox[
RowBox[{\(backward[x_?NumberQ]\), ":=",
RowBox[{"y", "/.",
RowBox[{"First", "[",
RowBox[{"NDSolve", "[",
RowBox[{
RowBox[{"{",
RowBox[{"de", ",", \(y(6000) == x\), ",",
RowBox[{
RowBox[{
SuperscriptBox["y", "\[Prime]",
MultilineFunction->None], "(", "6000", ")"}], "==",
"0"}]}], "}"}], ",", "y", ",", \({t, 10, 6000}\)}],

Cell[TextData[{
"supplying a final value, ",
Cell[BoxData[
". "
}], "Text"],

Cell[TextData[{
"Shooting backwards towards the initial condition, we use ",
Cell[BoxData[
FormBox[
StyleBox["FindRoot",
" to determine the final value, ",
Cell[BoxData[
":"
}], "Text"],

Cell[CellGroupData[{

Cell[BoxData[
1, 2}]\)], "Input"],

Cell[BoxData[
}, Open  ]],

Cell[TextData[{
"and hence find the (",
Cell[BoxData[
FormBox[
StyleBox["InterpolatingFunction",
") solution over the specified domain:"
}], "Text"],

Cell[CellGroupData[{

Cell[BoxData[
\(TraditionalForm\`bsoln = backward[x\ 10\^\(-10\)] /. %\)], "Input"],

Cell[BoxData[
FormBox[
TagBox[
RowBox[{"InterpolatingFunction", "[",
RowBox[{
RowBox[{"(", "\[NoBreak]", GridBox[{
{"10.`", "6000.`"}
},
ColumnAlignments->{Decimal}], "\[NoBreak]", ")"}],
",", "\<\"<>\"\>"}], "]"}],
False,
}, Open  ]],

Cell[TextData[{
"Note that the final value is ",
StyleBox["very",
FontSlant->"Italic"],
" small. ",
"Here we plot this solution:"
}], "Text"],

Cell[BoxData[
Plot[bsoln[t], {t, 10, 6000}, PlotRange -> {0, 30},
PlotStyle -> Hue];\)\)], "Input"],

Cell["\<\
Importantly, we can also see how well our numerical solution \
satisifies the differential equation:\
\>", "Text"],

Cell[BoxData[
Evaluate[First[de] /. y -> bsoln], {t, 10, 6000},
PlotRange -> {\(-10\^\(-6\)\), 10\^\(-6\)}];\)\)], "Input"],

Cell["\<\
The largest error, not suprisingly, is near our starting condition.\
\
\>", "Text"],

Cell["\<\
Next, we use our backward solution to determine the initial \
slope,\
\>", "Text"],

Cell[CellGroupData[{

Cell[BoxData[

Cell[BoxData[
}, Open  ]],

Cell["and hence the forward solution:", "Text"],

Cell[CellGroupData[{

Cell[BoxData[
FormBox[
RowBox[{"fsoln", "=",
RowBox[{"forward", "[",
FormBox[\(\(bsoln\^\[Prime]\)\),

Cell[BoxData[
FormBox[
TagBox[
RowBox[{"InterpolatingFunction", "[",
RowBox[{
RowBox[{"(", "\[NoBreak]", GridBox[{
{"10.`", "6000.`"}
},
ColumnAlignments->{Decimal}], "\[NoBreak]", ")"}],
",", "\<\"<>\"\>"}], "]"}],
False,
}, Open  ]],

Cell[TextData[{
"Plotting this solution together with the backward solution we see that \
they agree over ",
Cell[BoxData[
" but then diverge:"
}], "Text"],

Cell[BoxData[
Plot[fsoln[t], {t, 10, 6000}, PlotRange -> {0, 30},
PlotStyle -> Hue[1\/3], Epilog -> First[bplot]];\)\)], "Input"],

Cell["\<\
It is clear that the solution depends very sensitively on the \
initial (or final) values.  \
\>", "Text"],

Cell["\<\
Checking how well our forward numerical solution satisifies the \
differential equation we see that it is not as good as the backward solution:\

\
\>", "Text"],

Cell[BoxData[
Evaluate[First[de] /. y -> fsoln], {t, 10, 6000},
PlotRange -> {\(-10\^\(-6\)\), 10\^\(-6\)}];\)\)], "Input"],

Cell["\<\
In other words, for the boundary values you have supplied, this \
problem is not well conditioned. \
\>", "Text"]
}
]

```

• Prev by Date: Re: Verifying PrimeQ for n >10^16
• Next by Date: magic pen