```
In [1]:
```using Plots,LaTeXStrings
default(markersize=3,linewidth=1.5)
using DataFrames
using DifferentialEquations
include("FNC.jl");

```
In [2]:
```f = (u,t) -> sin((t+u)^2);
tspan = (0.0,4.0);
u0 = -1.0;

`DifferentialEquations`

package offers solvers for a huge variety of problems, including ordinary IVPs. Its syntax is more general (and powerful) than the one we will use to develop some simple solvers of our own.

```
In [3]:
```ivp = ODEProblem((u,p,t)->f(u,t),u0,tspan)
sol = solve(ivp);

```
In [4]:
```plot(sol,label="",ylabel="u(t)",title="Solution of u'=sin((t+u)^2)")

```
Out[4]:
```

`DifferntialEquations`

acts like any callable function that can be evaluated at different values of $t$, or plotted.

```
In [5]:
```f = (u,t) -> sin((t+u)^2);
tspan = (0.0,4.0);
u0 = -1.0;
ivp = ODEProblem((u,p,t)->f(u,t),u0,tspan)
sol = solve(ivp);

```
In [6]:
```@show sol(1.0);
@show sol.(0:.25:4);
plot(sol,label="u(t)")

```
Out[6]:
```

The object holds some information about how the values and plot are produced:

```
In [7]:
``````
sol
```

```
Out[7]:
```

```
In [8]:
```scatter!(sol.t,sol.u,label="discrete values")

```
Out[8]:
```

The equation $u'=(u+t)^2$ gives us some trouble.

```
In [9]:
```f = (u,t) -> (t+u)^2
using DifferentialEquations
sol = solve( ODEProblem((u,p,t)->f(u,t),1.,(0.,1.)) );

```
```

```
In [10]:
```using Plots
plot(sol,label="",yscale=:log10)

```
Out[10]:
```

```
In [11]:
```u = [ t->exp(t)*u0 for u0 in [0.7,1,1.3] ]
plot(u,0,3,leg=:none,xlabel="t",ylabel="u(t)", title="Exponential divergence of solutions")

```
Out[11]:
```

But with $u'=-u$, solutions actually get closer together with time.

```
In [12]:
```u = [ t->exp(-t)*u0 for u0 in [0.7,1,1.3] ]
plot(u,0,3,leg=:none,xlabel="t",ylabel="u(t)", title="Exponential convergence of solutions")

```
Out[12]:
```

We consider the IVP $u'=\sin[(u+t)^2]$ over $0\le t \le 4$, with $u(0)=-1$.

```
In [13]:
```f = (u,t) -> sin((t+u)^2);
tspan = (0.0,4.0);
u0 = -1.0;

```
In [14]:
```t,u = FNC.eulerivp(f,tspan,u0,20);
plot(t,u,m=:o,label=L"n=20",
xlabel="t", ylabel="u(t)", title="Solution by Euler's method" )

```
Out[14]:
```

```
In [15]:
```t,u = FNC.eulerivp(f,tspan,u0,200)
plot!(t,u,label=L"n=200")

```
Out[15]:
```

`DifferentialEquations`

solver to construct an accurate solution.

```
In [16]:
```ivp = ODEProblem((u,p,t)->f(u,t),u0,tspan)
u_exact = solve(ivp,reltol=1e-14,abstol=1e-14);
plot!(u_exact,l=:dash,label="accurate")

```
Out[16]:
```

Now we can perform a convergence study.

```
In [17]:
```n = @. 50*2^(0:5)
err = zeros(size(n))
for (j,n) = enumerate(n)
t,u = FNC.eulerivp(f,tspan,u0,n)
err[j] = maximum( @. abs(u_exact(t)-u) )
end
DataFrame(n=n,error=err)

```
Out[17]:
```

```
In [18]:
```plot(n,err,m=:o,label="results",
xaxis=(:log10,"n"), yaxis=(:log10,"inf-norm error"), title="Convergence of Euler's method")
plot!(n,0.05*(n/n[1]).^(-1),l=:dash,label="1st order")

```
Out[18]:
```

```
In [19]:
```A = [ -2 5; -1 0 ]

```
Out[19]:
```

```
In [20]:
```u0 = [1,0]
t = LinRange(0,6,600) # times for plotting
u = zeros(length(t),length(u0))
for j=1:length(t)
ut = exp(t[j]*A)*u0
u[j,:] = ut'
end

```
In [21]:
```plot(t,u,label=[L"u_1" L"u_2"],xlabel="t",ylabel="u(t)")

```
Out[21]:
```

We encode the predator–prey equations via a function.

```
In [22]:
```function predprey(u,p,t)
alpha,beta = p; y,z = u; # rename for convenience
s = (y*z) / (1+beta*y) # appears in both equations
return [ y*(1-alpha*y) - s, -z + s ]
end

```
Out[22]:
```

Note that the function accepts three inputs, `u`

, `p`

, and `t`

, even though there is no explicit dependence on `t`

. The second input is used to pass parameters that don't change throughout a single instance of the problem.

To solve the IVP we must also provide the initial condition, which is a 2-vector here, and the interval for the independent variable.

```
In [23]:
```u0 = [1,0.01]
tspan = (0.,80.)
alpha = 0.1; beta = 0.25;
sol = solve( ODEProblem(predprey,u0,tspan,[alpha,beta]) );

```
In [24]:
```plot(sol,label=["prey" "predator"],title="Predator-prey solution")

```
Out[24]:
```

`sol.u`

value is a vector of vectors, which can make manipulating the values a bit tricky. Here we convert the solution values to a matrix with two columns (one for each component).

```
In [25]:
```@show size(sol.u);
@show (sol.t[20],sol.u[20]);
u = [ sol.u[i][j] for i=1:length(sol.t), j=1:2 ]
plot!(sol.t,u,m=(:0,3),label="")

```
Out[25]:
```

*phase plane*, i.e., with $u_1$ and $u_2$ along the axes and time as a parameterization of the curve.

```
In [26]:
```plot(sol,vars=(1,2),label="",
xlabel="y",ylabel="z",title="Predator-prey phase plane")

```
Out[26]:
```

We consider the IVP $u'=\sin[(u+t)^2]$ over $0\le t \le 4$, with $u(0)=-1$.

```
In [27]:
```f = (u,t) -> sin((t+u)^2);
tspan = (0.0,4.0);
u0 = -1.0;

We use a `DifferentialEquations`

solver to construct an accurate approximation to the exact solution.

```
In [28]:
```ivp = ODEProblem((u,p,t)->f(u,t),u0,tspan)
u_exact = solve(ivp,reltol=1e-14,abstol=1e-14);

Now we perform a convergence study of our two Runge--Kutta implementations.

```
In [29]:
```n = @. 50*2^(0:5)
err_IE2 = zeros(size(n))
err_RK4 = zeros(size(n))
for (j,n) = enumerate(n)
t,u = FNC.ie2(f,tspan,u0,n)
err_IE2[j] = maximum( @.abs(u_exact(t)-u) )
t,u = FNC.rk4(f,tspan,u0,n)
err_RK4[j] = maximum( @.abs(u_exact(t)-u) )
end

```
In [30]:
```plot([2n 4n],[err_IE2 err_RK4],m=:o,label=["IE2" "RK4"],
xaxis=(:log10,"f-evaluations"),yaxis=(:log10,"inf-norm error"),
title="Convergence of RK methods",leg=:bottomleft)
plot!(2n,0.01*(n/n[1]).^(-2),l=:dash,label="2nd order")
plot!(4n,1e-6*(n/n[1]).^(-4),l=:dash,label="4th order")

```
Out[30]:
```

The fourth-order variant is more efficient in this problem over a wide range of accuracy.

Let's run adaptive RK on $u'=e^{t-u\sin u}$.

```
In [31]:
```f = (u,t) -> exp(t-u*sin(u))
t,u = FNC.rk23(f,[0.,5.],0.0,1e-5)
scatter(t,u,label="",
xlabel="t",ylabel="u(t)",title="Adaptive IVP solution")

```
Out[31]:
```

```
In [32]:
```plot(t[1:end-1],diff(t),label="",
xaxis=("t"),yaxis=(:log10,"step size"),title=("Adaptive step sizes"))

```
Out[32]:
```

If we had to run with a uniform step size to get this accuracy, it would be

```
In [33]:
```h_min = minimum(diff(t))

```
Out[33]:
```

On the other hand, the average step size that was actually taken was

```
In [34]:
```h_avg = sum(diff(t))/(length(t)-1)

```
Out[34]:
```

```
In [35]:
```f = (u,t) -> (t+u)^2;
t,u = FNC.rk23(f,[0,1],1.0,1e-5);

```
```

```
In [36]:
```plot(t,u,m=(:o,1),label="",
xlabel="t",yaxis=(:log10,"u(t)"),title="Finite-time blowup")

```
Out[36]:
```