Mathematica 9 is now available
Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2000
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2000

[Date Index] [Thread Index] [Author Index]

Search the Archive

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



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