# Taylor integration example: $\dot x=x^2$

Here, we will integrate the initial-value problem (IVP) defined by

$$
\begin{align}
\dot x &= x^2 \\
x(0)&=x_0
\end{align}
$$

Given a real number $A>0$, the restriction of the function $f(x)=x^2$ over the interval $I_A=(-A,A)$ satisfies the Lipschitz condition $|f(x_1)-f(x_2)|=|x_1^2-x_2^2|=|x_1+x_2|\cdot|x_1-x_2| \leq 2A|x_1-x_2|$ for every $x_1,x_2 \in I_A$. Therefore, the Picard-Lindelöf theorem for ordinary differential equations guarantees there exists a $\delta>0$ such that a solution for this IVP exists and is unique for $t \in [0,\delta]$, for any $x_0 \in I_A$. Note that in this case, it is necessary to restrict the function to a bounded interval in order to obtain a Lipschitz condition, which in turn is necessary to fulfill the conditions of the Picard-Lindelöf theorem.

Now, for any initial condition $x_0$, $t_0\in\mathbb{R}$ the analytical solution for this problem is:

$$
x(t)=\frac{x_0}{1-x_0\cdot (t-t_0)}
$$

In particular, the analytical solution exhibits a divergence at $t^*=t_0+1/x_0$; i.e., the solution is guaranteed to exist only for $t \in [t_0,t_0+1/x_0)$ if $x_0>0$. Otherwise, if $x_0<0$, then the analytical solution "moves away" from the singularity, since in this case we have $t^*=t_0+1/x_0<t_0$ (i.e., if $x_0<0$, the singularity is in the "past"). How does a numerical integrator behave near this divergence, when $x_0>0$?

We will try three methods to integrate this problem:

+ Adaptive time-step, 4th-order, Runge-Kutta (`ODE.jl`)
+ Taylor method (`TaylorIntegration.jl`)
+ Adaptive time-step, Runge-Kutta-Fehlberg 7/8 method (`ODE.jl`)

As initial conditions to integrate this IVP, we choose $x_0=3$, $t_0=0$, since then the singularity will be at $t^*=1/3$, and this number is not exactly representable in a binary floating-point format. Thus, any constant time-step numerical integrator should break down when integrating up to $t_\mathrm{max}=1/3$, or beyond.

We start off by including the relevant packages:


```julia
using TaylorIntegration, ODE, Plots, LaTeXStrings
```

The ODE:


```julia
diffeq(t, x) = x.^2
```




    diffeq (generic function with 1 method)



## 1. Adaptive time-step, 4th order, Runge Kutta method

We select $x_0=3$, $t_0=0$. Then, the singularity is at $t=1/3$.


```julia
@time tRK, xRK = ode45(diffeq, 3.0, [0.0, 0.34]); #warmup lap
@time tRK, xRK = ode45(diffeq, 3.0, [0.0, 0.34]);
```

    Warning: dt < minstep.  Stopping.
      1.397936 seconds (1.51 M allocations: 64.582 MB, 1.80% gc time)
    Warning: dt < minstep.  Stopping.
      0.000970 seconds (11.54 k allocations: 465.828 KB)


Plot $x$ vs $t$ (log-log):


```julia
pyplot()
plot(log10(tRK[2:end]), log10(xRK[2:end]))
title!("x vs t (log-log)")
xlabel!(L"\log_{10}(t)")
ylabel!(L"\log_{10}(x(t))")
```








    sys:1: MatplotlibDeprecationWarning: The set_axis_bgcolor function was deprecated in version 2.0. Use set_facecolor instead.


What is the final state of the system?


```julia
tRK[end], xRK[end]
```




    (0.33333423758781194,7.125446356124256e17)



Does the integrator get past the singularity?


```julia
tRK[end]>1/3
```




true



The answer is yes! So the last value of the solution is meaningless:


```julia
xRK[end] #this value is meaningless
```




7.125446356124256e17



How many steps did the RK integrator perform?


```julia
length(xRK)-1
```




166



How does the numerical solution compare to the analytical solution? The analytical solution is:


```julia
exactsol(t, x0) = x0./(1.0-x0.*t) #analytical solution
```




    exactsol (generic function with 1 method)



The relative difference between the numerical and analytical solution, $\delta x$, is:


```julia
δxRK = (xRK-exactsol(tRK, 3.0))./exactsol(tRK, 3.0) #numerical error, relative to analytical solution
;
```

The $\delta x$ vs $t$ plot (semilog):


```julia
plot(tRK[2:end], log10(abs(δxRK[2:end])))
title!("Relative error (semi-log)")
xlabel!(L"\log_{10}(t)")
ylabel!(L"\log_{10}(\delta x(t))")
```








This plot means that the error of the numerical solution grows systematically; and at the end of the integration, the error in the numerical solution is


```julia
(xRK[end-1]-exactsol(tRK[end-1], 3.0))
```




5.433622064235371e17



## 2. Taylor method

Again, we select $x_0=3$, $t_0=0$. The order of the Taylor integration is $28$, and we set the absolute tolerance equal to $10^{-20}$; this value is used during each time-step in order to compute an adaptive step size. We set the maximum number of integration steps equal to the number of steps that the previous integrator did.


```julia
@time tT, xT = taylorinteg(diffeq, 3.0, 0.0, 0.34, 28, 1e-20, maxsteps=length(xRK)-1); #warmup lap
@time tT, xT = taylorinteg(diffeq, 3.0, 0.0, 0.34, 28, 1e-20, maxsteps=length(xRK)-1);
```

    [1m[31mWARNING: Maximum number of integration steps reached; exiting.[0m


      0.215803 seconds (160.86 k allocations: 8.295 MB)
      0.001468 seconds (19.65 k allocations: 2.171 MB)



```julia
tT[end], xT[end]
```




    (0.3333333329479479,2.5948055925168757e9)



How many steps did the Taylor integrator perform?


```julia
length(xT)-1
```




166



Below, we show the $x$ vs $t$ plot (log-log):


```julia
plot(log10(tT[2:end]), log10(xT[2:end]))
title!("x vs t (log-log)")
xlabel!(L"\log_{10}(t)")
ylabel!(L"\log_{10}(x(t))")
```








Does the integrator get past the singularity?


```julia
tT[end] > 1/3
```




false



The answer is no! Even if we increase the value of the `maxsteps` keyword in `taylorinteg`, it doesn't get past the singularity!

Now, the relative difference between the numerical and analytical solution, $\delta x$, is:


```julia
δxT = (xT.-exactsol(tT, 3.0))./exactsol(tT, 3.0);
```

The $\delta x$ vs $t$ plot (logscale):


```julia
plot(tT[6:end], log10(abs(δxT[6:end])))
title!("Relative error (semi-log)")
xlabel!(L"t")
ylabel!(L"\log_{10}(\delta x(t))")
```








We observe that, while the execution time is ~10 times longer wrt 4th-order RK, the numerical solution obtained by the Taylor integrator stays within $10^{-12}$ of the analytical solution, for a same number of steps.

Now, that happens if we use a higher order Runge Kutta method to integrate this problem?

## 3. Runge-Kutta-Fehlberg 7/8 method

Here we use the Runge-Kutta-Fehlberg 7/8 method, included in `ODE.jl`, to integrate the same problem as before.


```julia
@time t78, x78 = ode78(diffeq, 3.0, [0.0, 0.34]); #warmup lap
@time t78, x78 = ode78(diffeq, 3.0, [0.0, 0.34]);
```

    Warning: dt < minstep.  Stopping.
      0.379769 seconds (345.22 k allocations: 15.063 MB, 2.40% gc time)
    Warning: dt < minstep.  Stopping.
      0.001199 seconds (12.76 k allocations: 494.219 KB)


Plot $x$ vs $t$ (log-log):


```julia
plot(log10(t78[2:end]), log10(x78[2:end]))
title!("x vs t (log-log)")
xlabel!(L"\log_{10}(t)")
ylabel!(L"\log_{10}(x(t))")
```








What is the final state of the system?


```julia
t78[end], x78[end]
```




    (0.3333336190078445,1.6711546943686948e18)



Does the integrator get past the singularity?


```julia
t78[end]>1/3
```




true



The answer is yes! So the last value of the solution is meaningless:


```julia
x78[end] #this value is meaningless
```




1.6711546943686948e18



How many steps did the RK integrator perform?


```julia
length(x78)-1
```




91



The relative difference between the numerical and analytical solution, $\delta x$, is:


```julia
δx78 = (x78-exactsol(t78, 3.0))./exactsol(t78, 3.0) #error relative to analytical solution
;
```

The $\delta x$ vs $t$ plot (semilog):


```julia
plot(t78[2:end], log10(abs(δx78[2:end])))
title!("Relative error (semi-log)")
xlabel!(L"t")
ylabel!(L"\log_{10}(\delta x(t))")
```








This time, the RKF 7/8 integrator is "only" twice as fast as the Taylor integrator, but the error continues to be greater than the error from the latter by several orders of magnitude.

## 4. Adaptive 4th-order RK, stringer tolerance

As a last example, we will integrate once again our problem using a 4th-order adaptive RK integrator, but imposing a stringer tolerance:


```julia
@time tRK_, xRK_ = ode45(diffeq, 3.0, [0.0, 0.34], abstol=1e-8, reltol=1e-8 ); #warmup lap
@time tRK_, xRK_ = ode45(diffeq, 3.0, [0.0, 0.34], abstol=1e-8, reltol=1e-8 );
;
```

    Warning: dt < minstep.  Stopping.
      0.127591 seconds (383.26 k allocations: 16.016 MB, 6.19% gc time)
    Warning: dt < minstep.  Stopping.
      0.029600 seconds (335.10 k allocations: 13.837 MB, 31.90% gc time)


Now, the integrator takes 10 times longer to complete the integration than the Taylor method.

Does it get past the singularity?


```julia
tRK_[end] > 1/3
```




true



Yes! So, once again, the last value reported by the integrator is completely meaningless. But, has it attained a higher precision than the Taylor method? Well, let's calculate once again the numerical error relative to the analytical solution:


```julia
δxRK_ = (xRK_-exactsol(tRK_, 3.0))./exactsol(tRK_, 3.0);
```

And now, let's plot this relative error vs time:


```julia
plot(tRK_[2:end], log10(abs(δxRK_[2:end])))
title!("Relative error (semi-log)")
xlabel!(L"t")
ylabel!(L"\log_{10}(\delta x(t))")
ylims!(-20,20)
```








The numerical error has actually gotten worse! `TaylorIntegration.jl` is indeed a really competitive package to integrate ODEs.


```julia

```

89dd4067c46a31cbead4a4c6be600c6084df7458
ipynb
Jupyter Notebook
examples/x-dot-equals-x-squared.ipynb
SebastianM-C/TaylorIntegration.jl
f3575ee1caba43e21312062d960613ec2ccba325
MIT
2016-09-22T22:32:12.000Z
2022-03-23T13:35:18.000Z
examples/x-dot-equals-x-squared.ipynb
SebastianM-C/TaylorIntegration.jl
f3575ee1caba43e21312062d960613ec2ccba325
MIT
2016-09-21T05:43:08.000Z
2022-03-15T02:55:17.000Z
examples/x-dot-equals-x-squared.ipynb
SebastianM-C/TaylorIntegration.jl
f3575ee1caba43e21312062d960613ec2ccba325
MIT
2016-09-24T04:37:11.000Z
2022-03-25T13:48:07.000Z
Qwen/Qwen-72B
1. YES
2. YES