"To numerically compute solutions of ordinary differential equations, we will use the package `OrdinaryDiffEq.jl`. This package is part of a larger package, `DifferentialEquations.jl`.\n",
"\n",
"For the documentation of these Julia packages on numerical integration, see \n",
"for example, on a step $[t_0, t_{1}]$, two approximate solutions at time $t_{1}$ are computed using two different schemes: $x_{1}$ and $\\hat{x}_{1}$, with local errors $O(h^p)$ and $O(h^{p+1})$, respectively (the global errors are $O(h^{p-1})$ and $O(h^p)$). This allows the local error of the less accurate scheme to be estimated using the difference $x_{1}-\\hat{x}_{1}$. The step size can then be adjusted so that the norm of this difference is below a given tolerance `Tol`, or such that \n",
"In practice, absolute and relative tolerances are considered component-wise, and we define: $\\mathrm{sc}_i = \\mathrm{abstol}_i + \\mathrm{reltol}_i \\times \\max(|x_{0i}|,|x_{1i}|)$. The step size is then calculated to satisfy \n",
"We focus here on the [simple pendulum](https://en.wikipedia.org/wiki/Pendulum_(mathematics)). The physical principles of classical mechanics give the governing equation for the motion:\n",
"where $\\ddot{\\alpha}(t)$ denotes the second derivative of the angle $\\alpha$ with respect to time $t$.\n",
"\n",
"- Using the state variable $ x = (x_1, x_2) = (\\alpha, \\dot{\\alpha}) $, write the function $ f $ to represent the differential equation in the form $ \\dot{x}(t) = f(t, x(t)) $.\n",
"- Implement this function below and execute the code. Based on your observations, do you see an oscillation or a rotation of the pendulum? \n",
"- Replace the initial angular velocity $ \\dot{\\alpha}(0) $ with 2. Comments?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\"\"\"\n",
" RHS of the IVP\n",
" x : state vector\n",
" λ : parameter vector\n",
" t : time variable. Here, time does not explicitly appear; the system is autonomous.\n",
"\"\"\"\n",
"function pendulum(x, λ, t)\n",
"\n",
" xdot = similar(x)\n",
" g, l, k, m = λ\n",
" \n",
" # TO COMPLETE/MODIFY\n",
" xdot[:] = zeros(size(x))\n",
" \n",
" return xdot \n",
" \n",
"end\n",
"\n",
"#\n",
"g = 9.81\n",
"l = 10\n",
"k = 0\n",
"m = 1\n",
"λ = [g, l, k, m] # constant parameters\n",
"\n",
"theta0 = pi/3\n",
"x0 = [theta0, 0] # initial state\n",
"\n",
"t0 = 0\n",
"tf = 3*pi*sqrt(l/g)*(1 + theta0^2/16 + theta0^4/3072) # 2*approximation of the period\n",
"tspan = (t0, tf) # initial and final instants\n",
"\n",
"prob = ODEProblem(pendulum, x0, tspan, λ, reltol = 1.e-8, abstol = 1.e-8) # defining the problem in Julia\n",
"sol = solve(prob) # numerical integration\n",
"\n",
"plot(sol, label = [\"x₁(t)\" \"x₂(t)\"]) # plotting the solution"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"**Phase Diagram.**\n",
"\n",
"- Run the code below and interpret the graph: where are the oscillations, rotations, and stable and unstable equilibrium points?\n",
"- Consider the case where $k=0.15$ (introducing friction). What happens?"
"# Equilibrium point and the Jacobian matrix at this point\n",
"xe = [0, 0]\n",
"A = dvdp(xe)\n",
"\n",
"# Eigenvalues of the Jacobian matrix\n",
"eigvals(A)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"<div style=\"width:95%;\n",
" margin:10px;\n",
" padding:8px;\n",
" color:white;\n",
" background-color: rgb(46, 109, 4);\n",
" border-radius:10px;\n",
" font-weight:bold;\n",
" font-size:1.5em;\n",
" text-align:center;\">\n",
"Predator-Prey Model (Lotka-Volterra Model)\n",
"</div>\n",
"\n",
"The differential equation considered is given by the [Lotka-Volterra predation equations](https://en.wikipedia.org/wiki/Lotka%E2%80%93Volterra_equations)\n",
"- $x(t)$ is the population of prey as a function of time;\n",
"- $y(t)$ is the population of predators as a function of time;\n",
"\n",
"The following parameters characterize the interactions between the two species:\n",
"\n",
"- $\\alpha$, intrinsic reproduction rate of prey (constant, independent of predator population);\n",
"- $\\beta$, mortality rate of prey due to encounters with predators;\n",
"- $\\delta$, reproduction rate of predators based on prey encountered and consumed;\n",
"- $\\gamma$, intrinsic mortality rate of predators (constant, independent of prey population);\n",
"\n",
"**Questions.**\n",
"\n",
"- The point $(0, 0)$ is clearly a stable equilibrium point. There is another equilibrium point, which one?\n",
"- Display the other equilibrium point along with some trajectories of the system in the phase plane, within the limits of the plot below. What do you observe about the trajectories?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# ------------------------------\n",
"# BEGIN TO COMPLETE/MODIFY THE RIGHT-HAND SIDE OF THE IVP\n",
"# ...\n",
"\n",
"# END TO COMPLETE/MODIFY THE RIGHT-HAND SIDE OF THE IVP\n",
"- By taking the state variable $x=(x_1,x_2)=(y, \\dot{y})$, write the function $f$ that allows you to express the differential equation in the form $\\dot{x}(t) = f(t,x(t))$.\n",
"- Code this function below and run the code with $\\omega_0 = \\omega/2$: are the solutions bounded?\n",
"- Replace the angular frequency to have $\\omega_0 = \\omega$. Comments.\n",
"\n",
"Note: See the Wikipedia page on the [resonance phenomenon](https://en.wikipedia.org/wiki/Resonance) for more information."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\"\"\"\n",
" Right-hand side of the IVP\n",
" x : state vector\n",
" λ : parameter vector\n",
" t : time variable. Here, time explicitly intervenes, the system is non-autonomous.\n",
To numerically compute solutions of ordinary differential equations, we will use the package `OrdinaryDiffEq.jl`. This package is part of a larger package, `DifferentialEquations.jl`.
For the documentation of these Julia packages on numerical integration, see
for example, on a step $[t_0, t_{1}]$, two approximate solutions at time $t_{1}$ are computed using two different schemes: $x_{1}$ and $\hat{x}_{1}$, with local errors $O(h^p)$ and $O(h^{p+1})$, respectively (the global errors are $O(h^{p-1})$ and $O(h^p)$). This allows the local error of the less accurate scheme to be estimated using the difference $x_{1}-\hat{x}_{1}$. The step size can then be adjusted so that the norm of this difference is below a given tolerance `Tol`, or such that
$$\frac{\|x_{1}-\hat{x}_{1}\|}{\mathrm{Tol}}<1.$$
In practice, absolute and relative tolerances are considered component-wise, and we define: $\mathrm{sc}_i = \mathrm{abstol}_i + \mathrm{reltol}_i \times \max(|x_{0i}|,|x_{1i}|)$. The step size is then calculated to satisfy
We will solve the initial value problem $\dot{x}(t) = 1 + x^p(t)$ with $p \ge 1$ and $x(0)=0$.
- Implement the function $f$ that reformulates the differential equation as $\dot{x}(t)=f(t,x(t))$.
- Determine the smallest value of $p \ge 1$ for which the solution blows up in finite time (with a precision of $0.1$) less than the value of 3.
%% Cell type:code id: tags:
``` julia
"""
Right-hand side of the IVP
x : state vector
p : scalar parameter
t : time variable. Here, time does not appear explicitly, so the system is autonomous.
"""
function f(x,p,t)
return0.0
end
#
x0=0.0
t0=0.0
tf=3
tspan=(t0,tf)# initial and terminal times
#
xspan=0:0.01:2
#
pf=plot(title="xᵖ")
px=plot(title="x(t, p)")
#
p=1
prob=ODEProblem(f,x0,tspan,p,reltol=1.e-8,abstol=1.e-8)# problem definition in Julia
sol=solve(prob)# numerical integration
pf=plot!(pf,xspan,x->(f(x,p,t0)-1),label="p=$p")# plot the right-hand side
px=plot!(px,sol,label="p=$p")# plot the solution
#
p=1# TO MODIFY
prob=ODEProblem(f,x0,tspan,p,reltol=1.e-8,abstol=1.e-8)# problem definition in Julia
sol=solve(prob)# numerical integration
pf=plot!(pf,xspan,x->(f(x,p,t0)-1),label="p=$p")# plot the right-hand side
px=plot!(px,sol,label="p=$p")# plot the solution
#
p=2
prob=ODEProblem(f,x0,tspan,p,reltol=1.e-8,abstol=1.e-8)# problem definition in Julia
sol=solve(prob)# numerical integration
pf=plot!(pf,xspan,x->(f(x,p,t0)-1),label="p=$p")# plotting the second member
px=plot!(px,sol,label="p=$p")# plotting the solution
#
plot!(pf,xlims=(0,2),ylims=(0,4))
plot!(px,xlims=(t0,tf),ylims=(0,100))
plot(pf,px,size=(650,350))
```
%% Cell type:markdown id: tags:
<div style="width:95%;
margin:10px;
padding:8px;
color:white;
background-color: rgb(46, 109, 4);
border-radius:10px;
font-weight:bold;
font-size:1.5em;
text-align:center;">
Simple Pendulum
</div>
We focus here on the [simple pendulum](https://en.wikipedia.org/wiki/Pendulum_(mathematics)). The physical principles of classical mechanics give the governing equation for the motion:
where $\ddot{\alpha}(t)$ denotes the second derivative of the angle $\alpha$ with respect to time $t$.
- Using the state variable $ x = (x_1, x_2) = (\alpha, \dot{\alpha}) $, write the function $ f $ to represent the differential equation in the form $ \dot{x}(t) = f(t, x(t)) $.
- Implement this function below and execute the code. Based on your observations, do you see an oscillation or a rotation of the pendulum?
- Replace the initial angular velocity $ \dot{\alpha}(0) $ with 2. Comments?
%% Cell type:code id: tags:
``` julia
"""
RHS of the IVP
x : state vector
λ : parameter vector
t : time variable. Here, time does not explicitly appear; the system is autonomous.
"""
function pendulum(x,λ,t)
xdot=similar(x)
g,l,k,m=λ
# TO COMPLETE/MODIFY
xdot[:]=zeros(size(x))
returnxdot
end
#
g=9.81
l=10
k=0
m=1
λ=[g,l,k,m]# constant parameters
theta0=pi/3
x0=[theta0,0]# initial state
t0=0
tf=3*pi*sqrt(l/g)*(1+theta0^2/16+theta0^4/3072)# 2*approximation of the period
tspan=(t0,tf)# initial and final instants
prob=ODEProblem(pendulum,x0,tspan,λ,reltol=1.e-8,abstol=1.e-8)# defining the problem in Julia
sol=solve(prob)# numerical integration
plot(sol,label=["x₁(t)""x₂(t)"])# plotting the solution
```
%% Cell type:markdown id: tags:
**Phase Diagram.**
- Run the code below and interpret the graph: where are the oscillations, rotations, and stable and unstable equilibrium points?
- Consider the case where $k=0.15$ (introducing friction). What happens?
%% Cell type:code id: tags:
``` julia
#
g=9.81
l=1.5
k=0
m=1
λ=[g,l,k,m]# constant parameters
plt=plot(xlabel="x₁",ylabel="x₂",legend=false)# initialize the plot
fortheta0in0:(2*pi)/10:2*pi
theta0_princ=theta0
tf=3*pi*sqrt(l/g)*(1+theta0_princ^2/16+theta0_princ^4/3072)# 2*approximation of the period
Below, we compute the eigenvalues of the Jacobian matrix of the right-hand side of the IVP at the equilibrium point $ x_e = (0, 0) $. Comments.
%% Cell type:code id: tags:
``` julia
# Jacobian of the vdp function
dvdp(x)=ForwardDiff.jacobian(x->vdp(x,[],0),x)
# Equilibrium point and the Jacobian matrix at this point
xe=[0,0]
A=dvdp(xe)
# Eigenvalues of the Jacobian matrix
eigvals(A)
```
%% Cell type:markdown id: tags:
<div style="width:95%;
margin:10px;
padding:8px;
color:white;
background-color: rgb(46, 109, 4);
border-radius:10px;
font-weight:bold;
font-size:1.5em;
text-align:center;">
Predator-Prey Model (Lotka-Volterra Model)
</div>
The differential equation considered is given by the [Lotka-Volterra predation equations](https://en.wikipedia.org/wiki/Lotka%E2%80%93Volterra_equations)
- $x(t)$ is the population of prey as a function of time;
- $y(t)$ is the population of predators as a function of time;
The following parameters characterize the interactions between the two species:
- $\alpha$, intrinsic reproduction rate of prey (constant, independent of predator population);
- $\beta$, mortality rate of prey due to encounters with predators;
- $\delta$, reproduction rate of predators based on prey encountered and consumed;
- $\gamma$, intrinsic mortality rate of predators (constant, independent of prey population);
**Questions.**
- The point $(0, 0)$ is clearly a stable equilibrium point. There is another equilibrium point, which one?
- Display the other equilibrium point along with some trajectories of the system in the phase plane, within the limits of the plot below. What do you observe about the trajectories?
%% Cell type:code id: tags:
``` julia
# ------------------------------
# BEGIN TO COMPLETE/MODIFY THE RIGHT-HAND SIDE OF THE IVP
# ...
# END TO COMPLETE/MODIFY THE RIGHT-HAND SIDE OF THE IVP
# ------------------------------
# Parameters
α=2/3
β=4/3
γ=1
δ=1
λ=[α,β,γ,δ]
t0=0
tf=20
tspan=(t0,tf)
plt=plot(xlabel="x₁",ylabel="x₂",legend=false)
# ------------------------------
# BEGIN TO COMPLETE/MODIFY THE DISPLAY
# ...
# END TO COMPLETE/MODIFY THE DISPLAY
# ------------------------------
plot!(xlims=(0,4),ylims=(0,2.5),size=(650,350))
```
%% Cell type:markdown id: tags:
<div style="width:95%;
margin:10px;
padding:8px;
color:white;
background-color: rgb(46, 109, 4);
border-radius:10px;
font-weight:bold;
font-size:1.5em;
text-align:center;">
Resonance Phenomenon
</div>
Let $\omega$ and $\omega_0$ be two strictly positive real numbers. Consider the differential equation
- By taking the state variable $x=(x_1,x_2)=(y, \dot{y})$, write the function $f$ that allows you to express the differential equation in the form $\dot{x}(t) = f(t,x(t))$.
- Code this function below and run the code with $\omega_0 = \omega/2$: are the solutions bounded?
- Replace the angular frequency to have $\omega_0 = \omega$. Comments.
Note: See the Wikipedia page on the [resonance phenomenon](https://en.wikipedia.org/wiki/Resonance) for more information.
%% Cell type:code id: tags:
``` julia
"""
Right-hand side of the IVP
x : state vector
λ : parameter vector
t : time variable. Here, time explicitly intervenes, the system is non-autonomous.
