NDSolve with Explicit, Implicit Euler and Trapezoidal method

Question

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. 

Answer 1

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", 
  "Coefficients" -> "ImplicitRungeKuttaRadauIIACoefficients", 
  "DifferenceOrder" -> 1, "ImplicitSolver" -> {"FixedPoint", 
  AccuracyGoal -> MachinePrecision, PrecisionGoal -> MachinePrecision, 
  "IterationSafetyFactor" -> 1}}};

With this definition, your (corrected)

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", 
  "Coefficients" -> "ImplicitRungeKuttaRadauIIACoefficients", 
  "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.

Answer 2

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.