{
"nbformat": 4,
"nbformat_minor": 0,
"metadata": {
"colab": {
"provenance": []
},
"kernelspec": {
"name": "python3",
"display_name": "Python 3"
},
"language_info": {
"name": "python"
}
},
"cells": [
{
"cell_type": "markdown",
"source": [
"#**Lab 5**\n",
"\n",
"**Name: (your name)**\n",
"\n"
],
"metadata": {
"id": "6VpZ01tyznHo"
}
},
{
"cell_type": "markdown",
"source": [
"Some remarks before you start:\n",
"* Import the Python packages below each time you start a session.\n",
"* Before submitting your report, remove all ouputs by clicking on **Edit🠦Clear all ouputs**. Then download your report by clicking on **File🠦Download🠦.ipynb**. This is the file you will submit on Canvas.\n",
"\n",
"In this lab, you will practice:\n",
"* solving ordinary differential equations (ODE) of first and second order using the finite difference method\n",
"* solving and simulating the motion of a pendulum\n",
"\n",
"\\begin{array}{|c| c|}\n",
" \\hline\n",
" \\text{Problems} & \\text{Points}\\\\\n",
" \\hline\\hline\n",
" 1 & 8 \\\\\n",
" \\hline\n",
" 2, 3 & 5 \\\\\n",
" \\hline\n",
" 4, 7 & 6 \\\\\n",
" \\hline\n",
" 5 & 4 \\\\\n",
" \\hline\n",
" 6 & 3 \\\\\n",
" \\hline\n",
" \\text{Readability of your report} & 3 \\\\\n",
" \\hline\n",
" \\text{Total: 7} & \\text{Total: 40} \\\\\n",
" \\hline\n",
"\\end{array}\n",
"\n"
],
"metadata": {
"id": "3cVoL2qvzufn"
}
},
{
"cell_type": "markdown",
"source": [
"##**I. Import necessary Python packages**\n",
"Execute the following code each time you start a session."
],
"metadata": {
"id": "bQagzUIHzzb_"
}
},
{
"cell_type": "code",
"source": [
"from matplotlib.pyplot import *\n",
"from matplotlib.animation import FuncAnimation as animation\n",
"from matplotlib import rc\n",
"from numpy import *\n",
"rc('animation', html='jshtml')"
],
"metadata": {
"id": "sGBvQoGSzoih"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"## **II. Finite difference method for first-order ODE**\n"
],
"metadata": {
"id": "yGFzT_PXz431"
}
},
{
"cell_type": "markdown",
"source": [
"Let us illustrate the method through an example. Consider the initial-value problem $$y'=y+x,\\ \\ y(-1)=1$$\n",
"You discretize the $x$-axis as $-1=x_0< x_1 < x_2 <...$ where $x_i-x_{i-1}=h$. The differential equation evaluated at $x=x_i$ takes the form\n",
"$$y'(x_i)=y(x_i)+x_i\\quad \\quad \\quad (1)$$\n",
"Write $y_i=y(x_i)$ and approximate the derivative $y'(x_i)$ using *forward difference*\n",
"$$y'(x_i)\\approx \\frac{y(x_i+h)-y(x_i)}{h}=\\frac{y_{i+1}-y_i}{h}$$\n",
"You get an approximation of the equation (1) as follows:\n",
"$$\\frac{y_{i+1}-y_i}{h}=y_i+x_i$$\n",
"which reduces to the *recursive formula*\n",
"$$y_{i+1}=y_i+h(y_i+x_i)\\quad \\quad \\quad (2)$$\n",
"The initial condition gives $y_0=-1$. Using the recursive formula $(2)$ repeatedly, you can compute $y_1, y_2, y_3, ...$\n",
"\n",
"Note: The method illustrated above uses forward difference approximation for the derivative. It is also called the *Euler's method*. In [Exercise 3](#backward), you will approximate the derivative using backward difference approximation."
],
"metadata": {
"id": "5ff0bdatZKUx"
}
},
{
"cell_type": "code",
"source": [
"# solve y' = y + x with initial condition y(x0)=y0\n",
"x0 = -1\n",
"y0 = 1\n",
"# with step size h\n",
"h = 0.1\n",
"# the number of steps\n",
"N = 20\n",
"# Array x = [x0,x1,...,xN]\n",
"x = linspace(x0,x0+N*h,N+1)\n",
"# Array y = [y0,y1,...,yN]\n",
"y = zeros(N+1)\n",
"y[0] = y0\n",
"for i in range(N):\n",
" y[i+1] = y[i] + h*(y[i] + x[i])\n",
"print(x)\n",
"print(y)\n",
"plot(x,y)\n",
"show()"
],
"metadata": {
"id": "kxIif3jP1xeP"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"You can solve the above initial value problem by hand using the Integrating Factor method. The exact solution is\n",
"$$y(x)=e^{x+1}-x-1$$\n",
"To compare the approximate solution with the true solution, you plot them together. In the code above, replace the line **plot(x,y)** with\n",
"\n",
"```\n",
"plot(x,y, label='approximate')\n",
"plot(x,e**(x+1)-x-1, label='exact')\n",
"legend()\n",
"```\n",
"and then re-run the code."
],
"metadata": {
"id": "sM6pvTUPipkz"
}
},
{
"cell_type": "markdown",
"source": [
"###**Exercise 1**\n",
"Consider the initial value problem\n",
"\n",
"$$y'+xy=x,\\quad y(-2)=3$$\n",
"\n",
"* Find the exact solution to this problem.\n",
"* Write the recursive formula according to the finite difference method.\n",
"* Find the approximate value of $y(0)$ according to Euler' method with step size $h=0.1$\n",
"* Graph on the same plot the exact solution together with the approximate solutions corresponding to $h=0.4,\\ 0.2,\\ 0.1,\\ 0.05$ on the interval $x\\in [-2,2]$. Comment on what you observe. *Note that you need to use different values of $N$ for different values of $h$ so that the interval is $[-2,2]$.*"
],
"metadata": {
"id": "DIQOwXFQHl-b"
}
},
{
"cell_type": "markdown",
"source": [
"###**Exercise 2**\n",
"Let $y=y(x)$ be a function satisfying\n",
"$$y'+y^2=x,\\ \\ \\ y(0)=0.7$$\n",
"* Write the recursive formula according to the finite difference method.\n",
"* Graph on the same plot the approximate solutions corresponding to $h=0.1$ and $h=0.01$ on the interval $x\\in[0,2]$.\n",
"* Find the maximum value of $y(x)$ on this interval."
],
"metadata": {
"id": "t-u_0lHaFAt9"
}
},
{
"cell_type": "markdown",
"source": [
"###**Exercise 3**\n",
"The purpose of this exercise is to help you see that forward difference approximation is not the only way to approximate the derivative. Consider the initial-value problem $$y'=y+x,\\ \\ y(-1)=1$$\n",
"You can approximate $y'(x_i)$ using backward difference approximation as follows:\n",
"$$y'(x_i)\\approx \\frac{y(x_i)-y(x_{i-1})}{h}=\\frac{y_{i}-y_{i-1}}{h}\\quad\\quad\\quad (3)$$\n",
"The equation $y'(x_i)=y(x_i)+x_i$ is now approximated by\n",
"$$\\frac{y_{i}-y_{i-1}}{h}=y_i+x_i$$\n",
"which reduces to a recursive formula:\n",
"$$y_i=\\frac{y_{i-1}+hx_i}{1-h}$$\n",
"With the initial condition $y_0=-1$, you can use this recursive formula to compute $y_1, y_2, y_3, ...$\n",
"* Rewrite the recursive formula $(4)$ with the index $i$ being shifted to $i+1$.\n",
"* Use step size $h=0.1$. Graph the approximate solution on the interval $x\\in[-1,1]$. *Hint: this should be a simple adjustment of the code provided above.*\n",
"* With step size $h=0.1$, graph on the same plot the exact solution together with the approximate solution using forward difference and the approximate solution using backward difference. Which method gives a better approximation?"
],
"metadata": {
"id": "bKyNLONGHTc_"
}
},
{
"cell_type": "markdown",
"source": [
"##**IV. Finite diffence method for second-order ODE**"
],
"metadata": {
"id": "5SE_QQYYISHu"
}
},
{
"cell_type": "markdown",
"source": [
"Let us illutrate the method through an example. Consider the initial value problem\n",
"$$y''+xy'+y=1,\\quad y(-1)=2,\\ y'(-1)=1$$\n",
"You discretize the $x$-axis as $-1=x_0< x_1 < x_2 <...$ where $x_i-x_{i-1}=h$. The differential equation evaluated at $x=x_i$ takes the form\n",
"$$y''(x_i)+x_iy'(x_i)+y(x_i)=1\\quad \\quad \\quad (4)$$\n",
"Write $y_i=y(x_i)$. Approximate the first derivative $y'(x_i)$ using *backward difference*\n",
"$$y'(x_i)\\approx \\frac{y(x_i)-y(x_i-h)}{h}=\\frac{y_{i}-y_{i-1}}{h}$$\n",
"and approximate the second derivative $y''(x_i)$ using *centered difference*\n",
"$$y''(x_i)\\approx \\frac{y(x_i+h)-2y(x_i)+y(x_{i}-h)}{h^2}=\\frac{y_{i+1}-2y_i+y_{i-1}}{h^2}$$\n",
"You get an approximation of the equation $(4)$ as follows:\n",
"$$\\frac{y_{i+1}-2y_i+y_{i-1}}{h^2}+x_i\\frac{y_i-y_{i-1}}{h}+y_i=1$$\n",
"which reduces to the *recursive formula*\n",
"$$y_{i+1}=(2-h^2-hx_i)y_i+(hx_i-1)y_{i-1}+h^2$$\n",
"If you know $y_0$ and $y_1$, you can use this recursive formula to compute $y_2, y_3, y_4, ...$ The initial condition $y(-1)=2$ implies $y_0=2$. To determine $y_1$, note that\n",
"\n",
"$$1=y'(-1)=y'(x_0)\\approx \\frac{y(x_1)-y(x_0)}{h}=\\frac{y_1-y_0}{h}.$$\n",
"\n",
"Thus, $y_1\\approx y_0+h=2+h$."
],
"metadata": {
"id": "bn9jxU__8wFq"
}
},
{
"cell_type": "code",
"source": [
"# solve y''+xy'+y=1 with initial condition y(x0)=a and y'(x0)=b\n",
"x0 = -1\n",
"a = 2\n",
"b = 1\n",
"# with step size h\n",
"h = 0.1\n",
"# the number of steps\n",
"N = 20\n",
"# Array x = [x0,x1,...,xN]\n",
"x = linspace(x0,x0+N*h,N+1)\n",
"# Array y = [y0,y1,...,yN]\n",
"y = zeros(N+1)\n",
"y[0] = a\n",
"y[1] = a+b*h\n",
"for i in range(1,N):\n",
" y[i+1] = (2-h**2-h*x[i])*y[i] + (h*x[i]-1)*y[i-1] + h**2\n",
"plot(x,y)\n",
"show()"
],
"metadata": {
"id": "rKoazOeiIctn"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"The above initial value problem has an exact solution\n",
"$$y(x)=1+\\exp\\left(\\frac{1-x^2}{2}\\right)$$\n",
"To compare the approximate solution with the true solution, you plot them together. In the code above, replace the line **plot(x,y)** with\n",
"\n",
"```\n",
"plot(x,y, label='approximate')\n",
"plot(x,1+e**((1-x**2)/2), label='exact')\n",
"legend()\n",
"```\n",
"and then re-run the code."
],
"metadata": {
"id": "0-QB6QpfTDVz"
}
},
{
"cell_type": "markdown",
"source": [
"###**Exercise 4**\n",
"The purpose of this exercise is to help you experiment with a different way to approximate the derivative. Consider the same initial value problem as above. Instead of using backward difference, you will use forward difference:\n",
"$$y'(x_i)\\approx \\frac{y_{i+1}-y_i}{h}$$\n",
"* Write the recursive formula using this new approximation of derivative.\n",
"* With step size $h=0.1$, find an approximation of $y(-0.2)$.\n",
"* With step size $h=0.1$, graph on the same plot the exact solution together with the approximate solution using forward difference and the approximate solution using backward difference. Which method gives a better approximation?"
],
"metadata": {
"id": "6QtOwQDPWZ1M"
}
},
{
"cell_type": "markdown",
"source": [
"###**Exercise 5**\n",
"A pendulum is a mechanical system consisting of a mass suspended from a fixed point by a string, allowing it to swing back and forth under the force of gravity. The position of the mass at time $t$ is completely determined by the angle $u=u(t)$ between the string and the vertical axis (figure below). The angle $u(t)$ is positive if the mass is on the right of the vertical axis, negative if on the left, and zero if on the vertical axis. If there is no friction, the angle satisfies the following second-order ODE (derived from Newton's Second Law):\n",
"\n",
"$$u''+\\frac{g}{L}\\sin u=0$$\n",
"\n",
"where $g\\approx 9.81\\,m/s^2$ is the grativational acceleration, and $L$ is the length of the string."
],
"metadata": {
"id": "KLJTj1VxflXW"
}
},
{
"cell_type": "markdown",
"source": [
""
],
"metadata": {
"id": "jLtk_xBdGGw4"
}
},
{
"cell_type": "markdown",
"source": [
"* Discretize the time as $0=t_0< t_1 < t_2 <...$ where $t_i-t_{i-1}=h$. Use the finite difference method to write a recursive formula for $u_i=u(t_i)$.\n",
"* Suppose that the length of the string is $L=1$ (meter), the initial angle is $u(0)=\\pi/3$ (radian), and the initial angular velocity is $u'(0)=0$ (rad/s). Use step size $h=0.1$ and number of steps $N=100$ to plot the graph of $u(t)$. *Hint: update the code provided above, replacing $y$ by $u$ and $x$ by $t$.*"
],
"metadata": {
"id": "1ISLOnVgGkCq"
}
},
{
"cell_type": "markdown",
"source": [
"###**Exercise 6**\n",
"With the $xy$ coordinate system indicated in the pendulum figure in [Exercise 5](#pendulum), you can see that the coordinates of the mass are $x=L\\sin u$ and $y=-L\\cos u$, where $u=u(t)$ is the angle between the string and the vertical axis.\n",
"\n",
"* Simulate the motion of the pendulum by running the code\n",
"```\n",
"fig = figure()\n",
"ax = fig.add_subplot(aspect='equal')\n",
"def animate(i):\n",
" ax.clear()\n",
" ax.set_xlim(-L*1.2, L*1.2)\n",
" ax.set_ylim(-L*1.2, L*1.2)\n",
" x = L*sin(u[i])\n",
" y = -L*cos(u[i])\n",
" ax.plot([0,x], [0,y], lw=3, c='black')\n",
" ax.add_patch(Circle([x,y], 0.08, fc='black'))\n",
"motion = animation(fig, animate, frames=N, repeat=True, interval=h*1000)\n",
"```\n",
"Then execute the following command in a new cell\n",
"```\n",
"motion\n",
"```\n",
"* To make the mass a little smaller (still circular shape) and change its color to red (while keeping the color of the string), how do you fix the above code?"
],
"metadata": {
"id": "_76xirr3Lk59"
}
},
{
"cell_type": "markdown",
"source": [
"###**Exercise 7**\n",
"In real life, all unforced pendulums will eventually stop due to some kind of friction (such as air resistance). Friction is a force that opposes the motion. The faster the pendulum moves, the stronger the resistance. Thus, you may think of the friction force as being proportional to the angular velocity, which is $u'$. The ODE is modified to\n",
"\n",
"$$u''+cu'+\\frac{g}{L}\\sin u=0$$\n",
"\n",
"where $c>0$ is a coefficient of friction.\n",
"* Discretize the time as $0=t_0< t_1 < t_2 <...$ where $t_i-t_{i-1}=h$. Use the finite difference method to write a recursive formula for $u_i=u(t_i)$.\n",
"* Suppose that the length of the string is $L=1$, the coefficient of friction is $c=0.5$, the initial angle is $u(0)=\\pi/3$, and the initial angular velocity is $u'(0)=0$. Use step size $h=0.1$ and number of steps $N=100$ to plot the graph of $u(t)$.\n",
"* Simulate the motion of the pendulum."
],
"metadata": {
"id": "lggls-NIQqyZ"
}
}
]
}