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[10] = 29
>
> y'[6000] = 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"}]}],
";"}], TraditionalForm]], "Input"],
Cell[TextData[{
"with boundary conditions ",
Cell[BoxData[
\(TraditionalForm\`y(10) == 29\)]],
" and ",
Cell[BoxData[
FormBox[
RowBox[{
RowBox[{
SuperscriptBox["y", "\[Prime]",
MultilineFunction->None], "(", "6000", ")"}], "==", "0"}],
TraditionalForm]]],
", where"
}], "Text"],
Cell[BoxData[
\(TraditionalForm\`p(
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[
\(TraditionalForm\`q(
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",
"Input"], TraditionalForm]]],
" cannot solve this boundary value problem directly. However it is possible
\
to solve such equations using ",
Cell[BoxData[
FormBox[
StyleBox["NDSolve",
"Input"], TraditionalForm]]],
"."
}], "Text"],
Cell[TextData[{
"Firstly, it is useful to graph the behaviour of ",
Cell[BoxData[
\(TraditionalForm\`p(t)\)]],
","
}], "Text"],
Cell[BoxData[
\(TraditionalForm\`\(Plot[p[t], {t, 10, 6000},
PlotRange -> All];\)\)], "Input"],
Cell[TextData[{
"and ",
Cell[BoxData[
\(TraditionalForm\`q(t)\)]],
":"
}], "Text"],
Cell[BoxData[
\(TraditionalForm\`\(Plot[q[t], {t, 10, 6000},
PlotRange -> All];\)\)], "Input"],
Cell[TextData[{
"We observe that ",
Cell[BoxData[
\(TraditionalForm\`p(t)\)]],
" is essentially flat for ",
Cell[BoxData[
\(TraditionalForm\`t \[GreaterTilde] 200\)]],
"."
}], "Text"],
Cell[TextData[{
"Since ",
Cell[BoxData[
FormBox[
StyleBox["NDSolve",
"Input"], TraditionalForm]]],
" 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}\)}],
"]"}], "]"}]}]}], TraditionalForm]], "Input"],
Cell[TextData[{
"supplying an initial slope, ",
Cell[BoxData[
\(TraditionalForm\`m\)]],
", 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}\)}],
"]"}], "]"}]}]}], TraditionalForm]], "Input"],
Cell[TextData[{
"supplying a final value, ",
Cell[BoxData[
\(TraditionalForm\`x\)]],
". "
}], "Text"],
Cell[TextData[{
"Shooting backwards towards the initial condition, we use ",
Cell[BoxData[
FormBox[
StyleBox["FindRoot",
"Input"], TraditionalForm]]],
" to determine the final value, ",
Cell[BoxData[
\(TraditionalForm\`x\)]],
":"
}], "Text"],
Cell[CellGroupData[{
Cell[BoxData[
\(TraditionalForm\`FindRoot[\(backward[x\ 10\^\(-10\)]\)[10] == 29, {x,
1, 2}]\)], "Input"],
Cell[BoxData[
\(TraditionalForm\`{x -> 1.4815130926261595`}\)], "Output"]
}, Open ]],
Cell[TextData[{
"and hence find the (",
Cell[BoxData[
FormBox[
StyleBox["InterpolatingFunction",
"Input"], TraditionalForm]]],
") 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,
Editable->False], TraditionalForm]], "Output"]
}, Open ]],
Cell[TextData[{
"Note that the final value is ",
StyleBox["very",
FontSlant->"Italic"],
" small. ",
"Here we plot this solution:"
}], "Text"],
Cell[BoxData[
\(TraditionalForm\`\(bplot =
Plot[bsoln[t], {t, 10, 6000}, PlotRange -> {0, 30},
PlotStyle -> Hue[0]];\)\)], "Input"],
Cell["\<\
Importantly, we can also see how well our numerical solution \
satisifies the differential equation:\
\>", "Text"],
Cell[BoxData[
\(TraditionalForm\`\(Plot[
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[
\(TraditionalForm\`\(bsoln\^\[Prime]\)[10]\)], "Input"],
Cell[BoxData[
\(TraditionalForm\`\(-0.10250604074425176`\)\)], "Output"]
}, Open ]],
Cell["and hence the forward solution:", "Text"],
Cell[CellGroupData[{
Cell[BoxData[
FormBox[
RowBox[{"fsoln", "=",
RowBox[{"forward", "[",
FormBox[\(\(bsoln\^\[Prime]\)[10]\),
"TraditionalForm"], "]"}]}], TraditionalForm]], "Input"],
Cell[BoxData[
FormBox[
TagBox[
RowBox[{"InterpolatingFunction", "[",
RowBox[{
RowBox[{"(", "\[NoBreak]", GridBox[{
{"10.`", "6000.`"}
},
ColumnAlignments->{Decimal}], "\[NoBreak]", ")"}],
",", "\<\"<>\"\>"}], "]"}],
False,
Editable->False], TraditionalForm]], "Output"]
}, Open ]],
Cell[TextData[{
"Plotting this solution together with the backward solution we see that \
they agree over ",
Cell[BoxData[
\(TraditionalForm\`\([0, 1000]\)\)]],
" but then diverge:"
}], "Text"],
Cell[BoxData[
\(TraditionalForm\`\(fplot =
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[
\(TraditionalForm\`\(Plot[
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"]
}
]