# NDSolve with Explicit, Implicit Euler and Trapezoidal method

I am using Mathematica 9.0 both on Linux and Windows and I would like to integrate the Van der Pol equation numerically using various techniques such as Explicit and Implicit Euler and Trapezoidal methods. I inspected the documentations and the examples, but I was unable to find how to implement the Implicit Euler and Trapezoidal methods although they are present on the help page. Can someone guide me through this? My notebook is as follows(excuse me on my inexperience in posting Mathematica code I used the method here):

s =NDSolve[{x''[t]==-x[t]+10*(1-x[t]^2)*x'[t], x[0]==2, x'[0]==0}, x, {t, 0, 20}, Method-> "ExplicitEuler", "StartingStepSize"-> 1/100]
{{x->InterpolatingFunction[{{0.,20.}},<>]}}
Plot[Evaluate[x[t]/.s],{t, 0, 20}, PlotRange->All]

p=NDSolve[{x''[t]==-x[t]+10*(1-x[t]^2)*x'[t], x[0]==2, x'[0]==0}, x, {t, 0, 20}, Method® "ImplicitEuler", "StartingStepSize"-> 1/100]
NDSolve::bdmtd: The value of the option Method -> ImplicitEuler is not a known built-in method, a symbol that could be a user-defined method, or a list with a name followed by method options.

• The first argument of NDSolve should be {x''[t] - x[t] + 10*(1 - x[t]^2)*x'[t] == 0, x[0] == 2, x'[0] == 0}. – bbgodfrey Jan 1 '15 at 20:24

Although Implicit Euler is described in the documentation, it may not be an implemented Method. In fact, the Wolfram discussion of the Lotka–Volterra Equation actually defines Backward or Implicit Euler, suggesting that it is not an implemented Method:

BackwardEuler = {"FixedStep", Method -> {"ImplicitRungeKutta",
"DifferenceOrder" -> 1, "ImplicitSolver" -> {"FixedPoint",
AccuracyGoal -> MachinePrecision, PrecisionGoal -> MachinePrecision,
"IterationSafetyFactor" -> 1}}};


p = x /. First@NDSolve[{x''[t] - x[t] + 10*(1 - x[t]^2)*x'[t] == 0, x[0] == 2,
x'[0] == 0}, x, {t, 0, 20}, Method -> BackwardEuler, StartingStepSize -> 1/100]


yields

although it does generate error messages, suggesting that the result may not be accurate for t > .14. I hope this helps.

Trapezoidal Method

Although the documentation is not very clear, I believe that the trapezoidal method can be implemented similarly to the backward Euler method:

Trapezoidal = {"FixedStep", Method -> {"ImplicitRungeKutta",
"Coefficients" -> "ImplicitRungeKuttaLobattoIIIACoefficients",
"DifferenceOrder" -> 1, "ImplicitSolver" -> {"FixedPoint",
AccuracyGoal -> MachinePrecision, PrecisionGoal -> MachinePrecision,
"IterationSafetyFactor" -> 1}}};
q = x /. First@NDSolve[{x''[t] - x[t] + 10*(1 - x[t]^2)*x'[t] == 0, x[0] == 2,
x'[0] == 0}, x, {t, 0, 20}, Method -> Trapezoidal, StartingStepSize -> 1/100];
Plot[q[t], {t, 0, 20}]


Here, error messages first occur for t > .18, suggesting that the answer may not be accurate for larger t. Moreover, the two solution curves in this Answer, while qualitatively similar, are not identical.

Update to Accommodate Revised Question

The revised equation (with my minor corrections) has a quite different solution:

s = x /. First@NDSolve[{x''[t] == -x[t] + 10*(1 - x[t]^2)*x'[t], x[0] == 2,
x'[0] == 0}, x, {t, 0, 20}, Method -> "ExplicitEuler"];
Plot[s[t], {t, 0, 20}]


Neither of the two implicit methods as previously configured can handle the abrupt change in slope of the solution t = 9.1. However, simply removing StartingStepSize-> 1/100 allows NDSolve sufficient flexibility to solve the equation.

BackwardEuler = {"FixedStep", Method -> {"ImplicitRungeKutta",
"DifferenceOrder" -> 1, "ImplicitSolver" -> {"FixedPoint",
AccuracyGoal -> MachinePrecision, PrecisionGoal -> MachinePrecision,
"IterationSafetyFactor" -> 1}}};
p = x /. First@NDSolve[{x''[t] == -x[t] + 10*(1 - x[t]^2)*x'[t], x[0] == 2,
x'[0] == 0}, x, {t, 0, 20}, Method -> BackwardEuler];
Plot[p[t], {t, 0, 20}]


and

Trapezoidal = {"FixedStep", Method -> {"ImplicitRungeKutta",
"Coefficients" -> "ImplicitRungeKuttaLobattoIIIACoefficients",
"DifferenceOrder" -> 1, "ImplicitSolver" -> {"FixedPoint",
AccuracyGoal -> MachinePrecision, PrecisionGoal -> MachinePrecision,
"IterationSafetyFactor" -> 1}}};
q = x /. First@NDSolve[{x''[t] == -x[t] + 10*(1 - x[t]^2)*x'[t], x[0] == 2,
x'[0] == 0}, x, {t, 0, 20}, Method -> Trapezoidal];
Plot[q[t], {t, 0, 20}]


both produce curves indistinguishable from that just above.

• Well thanks for the answer and how I can implement the regular trapezoidal rule that is thought in all introductory numerical method classes? It is also in the documentation you linked(Some Basic Methods). – Vesnog Jan 1 '15 at 21:48
• Please see addition to Answer. Note that other methods can be implemented in the same way. – bbgodfrey Jan 1 '15 at 22:16
• Thanks I will try it, I find a bit confusing since I only use Mathematica occasionally to check results of other software packages actually. The syntax can get pretty confusing, I would be happy to read the part of the document you refer to by the way to grasp the concepts. – Vesnog Jan 1 '15 at 22:18
• Just click on the links included in the Answer. Unfortunately, the documentation is rather sketchy. I use Mathematica a lot, but I still find parts of it confusing. That is the nature, I believe, of powerful, multi-purpose software. You might also wish to read [Wikipedia] (en.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods) for a discussion nonlinear Runga Kutta methods. – bbgodfrey Jan 1 '15 at 22:19
• Thanks once again I will have a look at them at a convenient time. – Vesnog Jan 1 '15 at 22:34

Using Picard iterations I get this series solution:

First we convert it to standard form $x'=f(x(t),t)$ and apply Picard:

Problem:

Solve $x^{\prime\prime}\left( t\right) +x\left( t\right) -10\left( 1-x\left( t\right) ^{2}\right) x^{\prime}\left( t\right) =0$ with $x\left( 0\right) =2,x^{\prime}\left( 0\right) =0$

Let $x_{1}=x,x_{2}=x^{\prime}\,,$ then $x_{1}^{\prime}=x_{2},x_{2}^{\prime }=-x_{1}+10\left( 1-x_{1}^{2}\right) x_{2}$, therefore

\begin{align*} \begin{pmatrix} x_{1}^{\prime}\\ x_{2}^{\prime} \end{pmatrix} & = \begin{pmatrix} x_{2}\\ -x_{1}+10\left( 1-x_{1}^{2}\right) x_{2} \end{pmatrix} \\ & =f\left( x\left( t\right) ,t\right) \end{align*}

Using Picard iterations

$$x^{k+1}=x^{0}+\int_{0}^{t}f\left( x^{k}\left( \tau\right) ,\tau\right) d\tau$$

For $k=0$, and using initial conditions $\begin{pmatrix} x_{1}\\ x_{2}% \end{pmatrix} ^{0}=$ $\begin{pmatrix} 2\\ 0 \end{pmatrix}$ we obtain

\begin{align*} \begin{pmatrix} x_{1}\\ x_{2}% \end{pmatrix} ^{1} & = \begin{pmatrix} x_{1}\\ x_{2} \end{pmatrix} ^{0}+\int_{0}^{t} \begin{pmatrix} x_{2}\\ -x_{1}+10\left( 1-x_{1}^{2}\right) x_{2} \end{pmatrix} ^{0}d\tau\\ & =% \begin{pmatrix} 2\\ 0 \end{pmatrix} +\int_{0}^{t} \begin{pmatrix} 0\\ -2+10\left( 1-4\right) 0 \end{pmatrix} d\tau= \begin{pmatrix} 2\\ 0 \end{pmatrix} + \begin{pmatrix} 0\\ -2t \end{pmatrix} = \begin{pmatrix} 2\\ -2t \end{pmatrix} \end{align*}

For $k=1$

\begin{align*} \begin{pmatrix} x_{1}\\ x_{2}% \end{pmatrix} ^{2} & =% \begin{pmatrix} x_{1}\\ x_{2}% \end{pmatrix} ^{0}+\int_{0}^{t} \begin{pmatrix} x_{2}\\ -x_{1}+10\left( 1-x_{1}^{2}\right) x_{2} \end{pmatrix} ^{1}d\tau\\ & = \begin{pmatrix} 2\\ 0 \end{pmatrix} +\int_{0}^{t} \begin{pmatrix} -2\tau\\ -2+10\left( 1-4\right) \left( -2\tau\right) \end{pmatrix} d\tau= \begin{pmatrix} 2\\ 0 \end{pmatrix} + \begin{pmatrix} -t^{2}\\ 2t\left( 15t-1\right) \end{pmatrix} =% \begin{pmatrix} 2-t^{2}\\ 30t^{2}-2t \end{pmatrix} \end{align*}

For $k=2$

\begin{align*} \begin{pmatrix} x_{1}\\ x_{2} \end{pmatrix} ^{3} & = \begin{pmatrix} x_{1}\\ x_{2} \end{pmatrix} ^{0}+\int_{0}^{t} \begin{pmatrix} x_{2}\\ -x_{1}+10\left( 1-x_{1}^{2}\right) x_{2} \end{pmatrix} ^{2}d\tau= \begin{pmatrix} 2\\ 0 \end{pmatrix} +\int_{0}^{t} \begin{pmatrix} 30\tau^{2}-2\tau\\ -\left( 2-\tau^{2}\right) +10\left( 1-\left( 2-\tau^{2}\right) ^{2}\right) \left( 30\tau^{2}-2\tau\right) \end{pmatrix} d\tau\\ & = \begin{pmatrix} 2-t^{2}+10t^{3}\\ -\frac{300}{7}t^{7}+\frac{10}{3}t^{6}+240t^{5}-20t^{4}-\frac{899}{3} t^{3}+30t^{2}-2t \end{pmatrix} \end{align*}

Hence after 3 iterations, solution is

$$x_{1}=2-t^{2}+10t^{3}$$

We can keep going, but lets use M to do few more:

m = 4;
s = {2, 0};
Do[
s = s /. t -> tao;
s = {2, 0} + Integrate[{s[[2]], -s[[1]] + 10 (1 - s[[1]]^2) s[[2]]}, {tao, 0, t}]
,
{k, 0, m}
];
Expand[s[[1]]];


 Plot[s[[1]], {t, 0, 20}]


I then verified this using Maple, which supports series solution, but I think it does not use Picard, but I get pretty close series to Maple's with Picard:

eq:= diff(x(t),t\$2)=-x(t)+10*(1-x(t)^2)*diff(x(t),t);
ic:= x(0)=2, D(x)(0)=0;
Order:=10;
sol:=dsolve({eq,ic},x(t),type='series',pt=0);


The only problem with picard, one needs to do more iterations than the few I did above to converge. Since it is doing symbolic integration, it will take much longer time to do. But numerical integration can be used instead to speed it up.

• I am grateful for your effort you put in this post but I want to only apply Explicit, Implicit Euler and the Trapezoidal rules in solving for the ODE I have presented in order check with the results I obtained in Matlab. By the way due to a copy paste error from the rtf file generated the ODE was not correct in my original post and edited it accordingly. It would be s =NDSolve[{x''[t]==-x[t]+10*(1-x[t]^2)*x'[t], x[0]==2, x'[0]==0}, x, {t, 0, 20}, Method-> "ExplicitEuler", "StartingStepSize"-> 1/100] {{x->InterpolatingFunction[{{0.,20.}},<>]}} – Vesnog Jan 2 '15 at 14:20