diff --git a/M2/TP-EDO.ipynb b/M2/TP-EDO.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..ab4c1045bd856f3e8dcea94cac65c2fe4f30d1b8 --- /dev/null +++ b/M2/TP-EDO.ipynb @@ -0,0 +1,853 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "<img src=\"https://gitlab.irit.fr/toc/etu-n7/julia/-/raw/main/assets/TSE_Logo_2019.png\" alt=\"TSE\" height=\"200\" style=\"margin-left:50px\"/>\n", + "\n", + "# Dynamical system and numerical integration\n", + "\n", + "- Date : 2024-2025\n", + "- Approximative duration : 1h15" + ] + }, + { + "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", + "Introduction\n", + "</div>\n", + "\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`.\n", + "\n", + "For the documentation of these Julia packages on numerical integration, see \n", + "[the DifferentialEquations.jl documentation](https://diffeq.sciml.ai/stable/).\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# load packages\n", + "using OrdinaryDiffEq\n", + "using ForwardDiff\n", + "using LinearAlgebra\n", + "using Plots" + ] + }, + { + "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", + "Example 1\n", + "</div>\n", + "\n", + "We will solve the initial value problem\n", + "$$(IVP)\\left\\{\\begin{array}{l}\n", + "\\dot{x}_1(t)=x_1(t)+x_2(t)+\\sin t\\\\\n", + "\\dot{x}_2(t)=-x_1(t)+3x_2(t)\\\\\n", + "x_1(0)=-9/25\\\\\n", + "x_2(0)=-4/25,\n", + "\\end{array}\\right.\n", + "$$\n", + "\n", + "- Write the function $f$ that expresses the differential equation in the form $\\dot{x}(t) = f(t, x(t))$.\n", + "- Verify that the solution to $(IVP)$ is\n", + "$$\\begin{align*}\n", + "x_1(t)&= (-1/25)(13\\sin t+9\\cos t)\\\\\n", + "x_2(t)&= (-1/25)(3\\sin t+4\\cos t).\n", + "\\end{align*}$$\n", + "- Implement the function $f$ and execute the code below. Any comments?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Note: \n", + "# An IVP (Initial Value Problem) is an ordinary differential equation \n", + "# (ODE in English, EDO in French) with an added initial condition.\n", + "\n", + "\"\"\"\n", + " Right-hand side of the IVP\n", + " x : state vector\n", + " λ : parameter vector\n", + " t : time variable. Here, time appears explicitly, so the system is non-autonomous.\n", + "\"\"\"\n", + "function example1(x, λ, t)\n", + "\n", + " xdot = similar(x)\n", + " \n", + " # TO COMPLETE/MODIFY\n", + " xdot[:] = zeros(size(x))\n", + "\n", + " return xdot\n", + "end\n", + "\n", + "λ = Float64[] # no parameters\n", + "\n", + "t0 = 0\n", + "tf = 10\n", + "tspan = (t0, tf) # integration interval\n", + "\n", + "x0 = [-9/25, -4/25] # initial condition\n", + "\n", + "prob = ODEProblem(example1, x0, tspan, λ) # definition of the IVP\n", + "sol = solve(prob) # solve the IVP\n", + "p1 = plot(sol, label = [\"x₁(t)\" \"x₂(t)\"], title = \"default options\") # plot the solution" + ] + }, + { + "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", + "Step Size Control\n", + "</div>\n", + "\n", + "Good numerical integrators include automatic step size adjustment. For \n", + "[Runge-Kutta methods](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods), \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", + "$$\\frac{\\|x_{1}-\\hat{x}_{1}\\|}{\\mathrm{Tol}}<1.$$\n", + "\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", + "\n", + "$$\\mathrm{err} = \\sqrt{\\frac{1}{n}\\sum_{i=1}^n\\left(\\frac{x_{1i}-\\hat{x}_{1i}}{\\mathrm{sc}_i}\\right)^2}<1.$$\n", + "\n", + "Reference: [Hairer, Nørsett, Wanner (2008) Solving Ordinary Differential Equations I Nonstiff Problems](https://link.springer.com/book/10.1007/978-3-540-78862-1).\n", + "\n", + "In `Julia`, default values can be modified: `reltol = 1e-3` and `abstol = 1e-6` \n", + "(when these values are scalars, the same quantities are used for all components).\n", + "\n", + "Perform numerical integration for: \n", + "- `reltol = 1e-3` and `abstol = 1e-6`\n", + "- `reltol = 1e-6` and `abstol = 1e-9`\n", + "- `reltol = 1e-10` and `abstol = 1e-15`\n", + "\n", + "and display the results including the solution." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# computation with default options\n", + "sol = solve(prob, reltol = 1e-3, abstol = 1e-6)\n", + "p1 = plot(sol, label = [\"x₁(t)\" \"x₂(t)\"], title = \"reltol=1e-3, abstol=1e-6\")\n", + "\n", + "# TO COMPLETE/MODIFY\n", + "#\n", + "p2 = plot()\n", + "\n", + "#\n", + "p3 = plot()\n", + "\n", + "# exact solution\n", + "T = t0:(tf-t0)/100:tf\n", + "sinT = sin.(T) # vectorized operation\n", + "cosT = cos.(T)\n", + "p4 = plot(T, [(-1/25)*(13*sinT + 9*cosT) (-1/25)*(3*sinT + 4*cosT)], label = [\"x₁(t)\" \"x₂(t)\"], title = \"Exact Solution\")\n", + "\n", + "# display\n", + "plot(p1, p2, p3, p4, size=(650, 350))" + ] + }, + { + "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", + "On the Finite-Time Blow-Up of Solutions\n", + "</div>\n", + "\n", + "We will solve the initial value problem $\\dot{x}(t) = 1 + x^p(t)$ with $p \\ge 1$ and $x(0)=0$. \n", + "\n", + "- Implement the function $f$ that reformulates the differential equation as $\\dot{x}(t)=f(t,x(t))$.\n", + "- 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", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\"\"\"\n", + " Right-hand side of the IVP\n", + " x : state vector\n", + " p : scalar parameter\n", + " t : time variable. Here, time does not appear explicitly, so the system is autonomous.\n", + "\"\"\"\n", + "function f(x, p, t)\n", + " return 0.0\n", + "end\n", + "\n", + "#\n", + "x0 = 0.0\n", + "t0 = 0.0\n", + "tf = 3\n", + "tspan = (t0, tf) # initial and terminal times\n", + "\n", + "#\n", + "xspan = 0:0.01:2\n", + "\n", + "#\n", + "pf = plot(title = \"xᵖ\")\n", + "px = plot(title = \"x(t, p)\")\n", + "\n", + "#\n", + "p = 1\n", + "prob = ODEProblem(f, x0, tspan, p, reltol = 1.e-8, abstol = 1.e-8) # problem definition in Julia\n", + "sol = solve(prob) # numerical integration\n", + "pf = plot!(pf, xspan, x -> (f(x, p, t0)-1), label = \"p=$p\") # plot the right-hand side\n", + "px = plot!(px, sol, label = \"p=$p\") # plot the solution\n", + "\n", + "#\n", + "p = 1 # TO MODIFY\n", + "prob = ODEProblem(f, x0, tspan, p, reltol = 1.e-8, abstol = 1.e-8) # problem definition in Julia\n", + "sol = solve(prob) # numerical integration\n", + "pf = plot!(pf, xspan, x -> (f(x, p, t0)-1), label = \"p=$p\") # plot the right-hand side\n", + "px = plot!(px, sol, label = \"p=$p\") # plot the solution\n", + "\n", + "#\n", + "p = 2\n", + "prob = ODEProblem(f, x0, tspan, p, reltol = 1.e-8, abstol = 1.e-8) # problem definition in Julia\n", + "sol = solve(prob) # numerical integration\n", + "pf = plot!(pf, xspan, x -> (f(x, p, t0)-1), label = \"p=$p\") # plotting the second member\n", + "px = plot!(px, sol, label = \"p=$p\") # plotting the solution\n", + "\n", + "#\n", + "plot!(pf, xlims=(0, 2), ylims=(0, 4))\n", + "plot!(px, xlims=(t0, tf), ylims=(0, 100))\n", + "plot(pf, px, size=(650, 350))" + ] + }, + { + "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", + "Simple Pendulum\n", + "</div>\n", + "\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", + "\n", + "$$ ml^2\\ddot{\\alpha}(t)+mlg\\sin(\\alpha(t)) +k\\dot{\\alpha}(t)=0, $$\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?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#\n", + "g = 9.81\n", + "l = 1.5\n", + "k = 0\n", + "m = 1\n", + "λ = [g, l, k, m] # constant parameters\n", + "\n", + "plt = plot(xlabel = \"x₁\", ylabel = \"x₂\", legend = false) # initialize the plot\n", + "\n", + "for theta0 in 0:(2*pi)/10:2*pi\n", + " theta0_princ = theta0\n", + " tf = 3*pi*sqrt(l/g)*(1 + theta0_princ^2/16 + theta0_princ^4/3072) # 2*approximation of the period\n", + " tspan = (0.0,tf)\n", + " x0 = [theta0 0]\n", + " prob = ODEProblem(pendulum, x0, tspan, λ, reltol = 1.e-8, abstol = 1.e-8)\n", + " sol = solve(prob)\n", + " plot!(plt, sol, idxs=(1,2), color=\"blue\") # lw = linewidth\n", + "end\n", + "\n", + "theta0 = pi-10*eps()\n", + "x0 = [theta0 0]\n", + "tf = 50 # problem for tf=50 1/4 of the period!\n", + "tspan = (0.0,tf)\n", + "prob = ODEProblem(pendulum, x0, tspan, λ, reltol = 1.e-8, abstol = 1.e-8)\n", + "sol = solve(prob)\n", + "plot!(plt, sol, idxs=(1,2), xlims = (-2*pi,4*pi), color=\"green\") # lw = linewidth\n", + "\n", + "theta0 = pi+10*eps()\n", + "x0 = [theta0 0]\n", + "tf = 50\n", + "tspan = (0.0,tf)\n", + "prob = ODEProblem(pendulum, x0, tspan, λ, reltol = 1.e-8, abstol = 1.e-8)\n", + "sol = solve(prob)\n", + "plot!(plt, sol, idxs=(1,2), xlims = (-2*pi,4*pi), color=\"green\") # lw = linewidth\n", + "\n", + "# circulation case\n", + "for thetapoint0 in 0:1.:4 \n", + " tf = 10\n", + " tspan = (0.,tf)\n", + " x0 = [-pi thetapoint0] # thetapoint0 > 0 so theta increases from -pi to ...\n", + " prob = ODEProblem(pendulum, x0, tspan, λ, reltol = 1.e-8, abstol = 1.e-8)\n", + " sol = solve(prob)\n", + " plot!(plt, sol, idxs=(1,2), color=\"red\") # lw = linewidth\n", + "end\n", + "for thetapoint0 in -4:1.:0\n", + " tf = 10\n", + " tspan = (0.,tf)\n", + " x0 = [3*pi thetapoint0] # thetapoint0 < 0 so theta decreases from 3pi to ...\n", + " prob = ODEProblem(pendulum, x0, tspan, λ, reltol = 1.e-8, abstol = 1.e-8)\n", + " sol = solve(prob)\n", + " plot!(plt, sol, idxs=(1,2), color=\"purple\") # lw = linewidth\n", + "end\n", + "plot!(plt, [-pi 0 pi 2*pi 3*pi], [0 0 0 0 0], seriestype=:scatter)\n", + "plot!(plt, xlims = (-pi,3*pi), size=(650, 350))" + ] + }, + { + "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", + "Van der Pol Model\n", + "</div>\n", + "\n", + "The differential equation considered is the [Van der Pol equation](https://en.wikipedia.org/wiki/Van_der_Pol_oscillator)\n", + "$$(IVP)\\left\\{\\begin{array}{l}\n", + "\\dot{x}_1(t)=x_2(t)\\\\\n", + "\\dot{x}_2(t)=(1-x_1^2(t))x_2(t)-x_1(t)\\\\\n", + "x_1(0)=2.00861986087484313650940188\\\\\n", + "x_2(0)=0\n", + "\\end{array}\\right.\n", + "$$\n", + "up to the time $t_f=T=6.6632868593231301896996820305$, where $T$ is the period of the solution.\n", + "\n", + "Code the right-hand side of the IVP below and execute the code. Comments." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\"\"\"\n", + " Right-hand side of the IVP: Van der Pol model\n", + " x : state vector\n", + " λ : parameter vector\n", + " t : time variable. Here, time does not explicitly intervene, the system is autonomous.\n", + "\"\"\"\n", + "function vdp(x, λ, t)\n", + "\n", + " xpoint = similar(x)\n", + " \n", + " # TO BE COMPLETED/MODIFIED\n", + " xpoint[:] = zeros(size(x))\n", + " \n", + " return xpoint\n", + " \n", + "end\n", + "\n", + "#\n", + "t0 = 0\n", + "tf = 6.6632868593231301896996820305\n", + "tspan = (t0, tf)\n", + "\n", + "#\n", + "plot(xlabel = \"x₁\", ylabel = \"x₂\", legend = false)\n", + "\n", + "# Inner trajectories\n", + "for x01 in -2:0.4:0\n", + " x0 = [x01,0]\n", + " prob = ODEProblem(vdp, x0, tspan, reltol = 1.e-8, abstol = 1.e-8)\n", + " sol = solve(prob)\n", + " plot!(sol, idxs=(1,2), color = \"blue\") \n", + "end\n", + "\n", + "# Outer trajectories\n", + "for x02 in 2.5:0.5:4\n", + " x0 = [0,x02]\n", + " prob = ODEProblem(vdp, x0, tspan, reltol = 1.e-8, abstol = 1.e-8)\n", + " sol = solve(prob)\n", + " plot!(sol, idxs=(1,2), color = \"green\")\n", + "end\n", + "\n", + "# Periodic limit cycle\n", + "x0 = [2.00861986087484313650940188,0]\n", + "prob = ODEProblem(vdp, x0, tspan, reltol = 1.e-8, abstol = 1.e-8)\n", + "sol = solve(prob)\n", + "plot!(sol, idxs=(1,2), color = \"red\", lw = 2.)\n", + "\n", + "#\n", + "plot!([0], [0], seriestype=:scatter) # equilibrium point\n", + "plot!(xlims = (-2.5, 2.5), ylims = (-3, 5), size=(650, 350))\n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Stability.**\n", + "\n", + "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", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Jacobian of the vdp function\n", + "dvdp(x) = ForwardDiff.jacobian(x -> vdp(x, [], 0), x)\n", + "\n", + "# 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", + "\n", + "$$\\left\\{\\begin{array}{l}\n", + "\\dot{x}_1(t)= \\phantom{-}x_1(t) ( \\alpha - \\beta x_2(t)) \\\\[0.5em]\n", + "\\dot{x}_2(t)= -x_2(t) ( \\gamma - \\delta x_1(t))\n", + "\\end{array}\\right.\n", + "$$\n", + "\n", + "where\n", + "\n", + "- $t$ is time;\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", + "# ------------------------------\n", + "\n", + "# Parameters\n", + "α = 2/3\n", + "β = 4/3\n", + "γ = 1\n", + "δ = 1\n", + "λ = [α, β, γ, δ]\n", + "\n", + "t0 = 0\n", + "tf = 20\n", + "tspan = (t0, tf)\n", + "\n", + "plt = plot(xlabel = \"x₁\", ylabel = \"x₂\", legend = false)\n", + "\n", + "# ------------------------------\n", + "# BEGIN TO COMPLETE/MODIFY THE DISPLAY\n", + "# ...\n", + "\n", + "# END TO COMPLETE/MODIFY THE DISPLAY\n", + "# ------------------------------\n", + "\n", + "plot!(xlims = (0, 4), ylims = (0, 2.5), size=(650, 350))" + ] + }, + { + "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", + "Resonance Phenomenon\n", + "</div>\n", + "\n", + "Let $\\omega$ and $\\omega_0$ be two strictly positive real numbers. Consider the differential equation\n", + "\n", + "$$\\ddot{y}(t) + \\omega_0^2 y(t) = \\cos(\\omega t).$$\n", + "\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", + "\"\"\"\n", + "function resonance(x, λ, t)\n", + "\n", + " xdot = similar(x)\n", + " \n", + " # TO BE COMPLETED/MODIFIED\n", + " xdot[:] = zeros(size(x))\n", + " \n", + " return xdot\n", + "\n", + "end\n", + "\n", + "# parameters\n", + "ω = 1\n", + "ω₀ = ω/2\n", + "λ = [ω, ω₀]\n", + "\n", + "# integration times\n", + "t0 = 0\n", + "tf = 100\n", + "tspan = (t0, tf)\n", + "\n", + "# initial condition\n", + "x0 = [0, 0]\n", + "\n", + "# define the problem and solve it\n", + "prob = ODEProblem(resonance, x0, tspan, λ, reltol = 1.e-8, abstol = 1.e-8)\n", + "sol = solve(prob)\n", + "\n", + "function makegif()\n", + "\n", + " #\n", + " plt = plot(xlabel = \"x₁\", ylabel = \"x₂\", zlabel = \"t\", legend = false)\n", + "\n", + " #\n", + " ω == ω₀ ? lim = 50 : lim = 5\n", + " plot!(xlims = (-lim, lim), ylims = (-lim, lim), zlims = (t0, tf), size=(850, 550))\n", + "\n", + " # get x1, x2, t from sol\n", + " x1 = [sol.u[i][1] for i ∈ 1:length(sol.u)]\n", + " x2 = [sol.u[i][2] for i ∈ 1:length(sol.u)]\n", + " t = sol.t\n", + "\n", + " # plot the trajectory step by step and make a gif\n", + " i = 1\n", + " @gif for j ∈ 2:2:length(sol.t)\n", + " # plot the trajectory part from i to j\n", + " plot!(plt, x1[i:j], x2[i:j], t[i:j], color=:blue)\n", + " # update i \n", + " i = j\n", + " end\n", + "\n", + "end\n", + "\n", + "makegif()" + ] + }, + { + "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", + "Lorenz Attractor\n", + "</div>\n", + "\n", + "The differential equation considered is the [Lorenz equation](https://en.wikipedia.org/wiki/Lorenz_system), given by\n", + "\n", + "$$\\left\\{\\begin{array}{l}\n", + "\\dot{x}_1(t)= \\phantom{-}\\sigma (x_2(t) - x_1(t)) \\\\[0.5em]\n", + "\\dot{x}_2(t)= \\rho\\, x_1(t) - x_2(t) - x_1(t)\\, x_3(t) \\\\[0.5em]\n", + "\\dot{x}_3(t)= x_1(t)\\, x_2(t) - \\beta x_3(t)\n", + "\\end{array}\\right.\n", + "$$\n", + "\n", + "where \n", + "\n", + "- $t$ is time;\n", + "- $x(t)$ is proportional to the intensity of the convection motion;\n", + "- $y(t)$ is proportional to the temperature difference between the rising and descending currents;\n", + "- $z(t)$ is proportional to the deviation of the vertical temperature profile from the reference linear value;\n", + "- $\\sigma$, $\\beta$, and $\\rho$ are positive real parameters;\n", + "- $\\sigma$ is the Prandtl number;\n", + "- $\\rho$ is the Rayleigh number.\n", + "\n", + "It is common to set $\\sigma = 10$, $\\beta = 8/3$, and $\\rho = 28$ to illustrate the phenomenon of chaos.\n", + "\n", + "Complete the code below and execute it. Observe the phenomenon of chaos. You can modify the initial condition." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\"\"\"\n", + " Second member of the IVP\n", + " x : state vector\n", + " λ : parameter vector\n", + " t : time variable. Here time does not explicitly intervene, the system is autonomous.\n", + "\"\"\"\n", + "function Lorentz(x, λ, t)\n", + "\n", + " xdot = similar(x)\n", + " σ, ρ, β = λ\n", + " \n", + " # TO COMPLETE/MODIFY\n", + " xdot[:] = zeros(size(x))\n", + " \n", + " return xdot\n", + " \n", + "end\n", + "\n", + "# parameters\n", + "σ = 10\n", + "ρ = 28\n", + "β = 8/3\n", + "λ = [σ, ρ, β]\n", + "\n", + "# integration times\n", + "t0 = 0\n", + "tf = 100\n", + "tspan = (t0, tf)\n", + "\n", + "# initial condition\n", + "x0 = [0, 1, 0]\n", + "\n", + "# define the problem and solve it\n", + "prob = ODEProblem(Lorentz, x0, tspan, λ, reltol = 1.e-8, abstol = 1.e-8)\n", + "sol = solve(prob)\n", + "\n", + "function makegif()\n", + "\n", + " #\n", + " plt = plot(xlabel = \"x₁\", ylabel = \"x₂\", zlabel = \"x₃\", legend = false)\n", + "\n", + " # \n", + " plot!(plt, background_color = :royalblue4)\n", + "\n", + " # update the color of axis of the figure to white\n", + " plot!(plt, axiscolor = :white, tickfontcolor = :white, bordercolor = :white, fg_color_guide = :white)\n", + "\n", + " # plot the three equilibrium points\n", + " plot!(plt, [0], [0], [0], seriestype=:scatter, color=:red)\n", + " plot!(plt, [ sqrt(β*(ρ-1))], [ sqrt(β*(ρ-1))], [ρ-1], seriestype=:scatter, color=:red)\n", + " plot!(plt, [-sqrt(β*(ρ-1))], [-sqrt(β*(ρ-1))], [ρ-1], seriestype=:scatter, color=:red)\n", + "\n", + " # get x1, x2, t from sol\n", + " x1 = [sol.u[i][1] for i ∈ 1:length(sol.u)]\n", + " x2 = [sol.u[i][2] for i ∈ 1:length(sol.u)]\n", + " x3 = [sol.u[i][3] for i ∈ 1:length(sol.u)]\n", + "\n", + " #\n", + " plot!(xlims = (-20, 20), ylims = (-25, 30), zlims = (0, 50), size=(850, 550))\n", + "\n", + " # \n", + " duration = 10\n", + " fps = 10\n", + " frames = fps*duration\n", + " step = length(sol.t) ÷ frames\n", + "\n", + " # return if there is nothing to compute\n", + " if (step == 0)\n", + " return;\n", + " end\n", + "\n", + " # plot the trajectory step by step and make a gif\n", + " i = 1\n", + " animation = @animate for j ∈ 2:step:length(sol.t)\n", + " # plot the trajectory part from i to j\n", + " plot!(plt, x1[i:j], x2[i:j], x3[i:j], color=:gold2)\n", + " # update i \n", + " i = j\n", + " end\n", + "\n", + " gif(animation, \"lorentz.gif\", fps=fps)\n", + "\n", + "end\n", + "\n", + "makegif()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.11.1", + "language": "julia", + "name": "julia-1.11" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.11.1" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}