Solving ODEs Using Python

Consider the initial value problem \(\frac{dx}{dt}=f(t,x), x(0)=x_0\). We will use the odeint solver from scipy library to solve this system.

Example 1: Newton's Law of Cooling

We will first solve the following initial value problem that models the Newton's law of cooling.

where \(x=x(t)\) is the temperature of the body at time \(t\), \(A\) is the temperature of its sorroundings, and K is a constant of proporionality which will depend on the initial temperature. For our computation purpose, we will take \(A=40\), initial temperature \(x_0=60\) degrees, and \(K=-0.03\). Then our initial value problem takes the form $$\frac{dx}{dt}=-0.03(x-40), \quad x(0)=60$$

Exercise 1

Assume that air resistance is proportional to the square of the instanteneous velocity. The velocity \(v\) of an object with mass \(m \) dropped from a given height is determined by $$ m \frac{dv}{dt}=mg-kv^2, \quad k>0, $$ where \( g\) is the acceleration due to gravity and k is the constant of proportionlality. Let \(v(0) =0, k=0.24, m=10 \) lb, and \( g=32 ft/s^2\).

1. Use the above method to approximate the velocity of the falling object after 6 sec.
2. Graph the numerical solution in the interval \( [0, 6] \).
3. Find the exact solution of this problem and draw the exact solution on the same frame.
4. Calculate the absolute error of the numerical solution at \( t=6\).

Example 2: Irregular Heartbeats and Lidocaine

The irregular heartbeat is treated clinically with lidocaine. The system of differential equations modelling this therapy is given as a system of two linear differential equations.

We can think of this system as $$\frac{d\textbf{y}}{dt}=\textbf{f}(x,\textbf{y}),$$ where \( \textbf{f}(x,\textbf{y})=(-0.09y_1+0.038y_2,0.066y_1-0.038y_2) \).

Exercise 2: Predator-Prey Model

As our next step, we will solve a system of two nonlinear Predator-Prey model, also known as the Lotka-Volterra system. This system is used to model the dynamics of biological system between two competing species, predator and prey. $$ y_1'=\alpha y_1-\beta y_1y_2$$ $$ y_2'=-\gamma y_2+\delta y_1y_2$$ with the initial conditions $$ y_1(0)=y_10, \quad y_2(0)=y_20, $$ where the scalar \(x \) is the time and the vector \(\textbf{y}=(y_1,y_2) \). The quantities \(y_1(x) \) and \(y_2(x) \) represent the population density of the predator and prey at time \(x\). The parameters, \( \alpha, \beta, \gamma \), and \( \delta \) represent the growth and death rates of predator and prey. All parameters are positive and real.

We can think of this system as $$\frac{d\textbf{y}}{dx}=\textbf{f}(x,\textbf{y})$$ where \( \textbf{f}(x,\textbf{y})=(\alpha y_1-\beta y_1y_2,-\gamma y_2+\delta y_1y_2) \).

Complete the following steps.

1. Use the codes with some modification to numerically solve the Lotka-Volterra system in the time interval [0, 18]. Use \( \alpha=2, \beta=1, \delta=3, \) and \( \gamma=9 \).
2. Save the data in a file as a 18x3 matrix or table with columns representing the value of \( t \), \( y_1 \) and \( y_2 \) respectively.
3. Draw the graphs of \( y_1 \) and \( y_2 \) versus t in the same frame.
4. Draw the graph of predator versus pray.
5. Identify the
6. Repeat the above steps using different values of the parameters and explain the pattern in the dynamics between the two species.

Example 3: Lorenz System

We will now solve the Lorenz system: $$\frac{dx}{dt}=\sigma(y-x)$$ $$\frac{dy}{dt}=x(\rho-z)-y $$ $$ \frac{dz}{dt}=xy-\beta z$$ The system is related to the properties of a two-dimensional fluid layer. The fluid is uniformly heated from below and cooled from above. The scalars \(x, y \) and \( z\) are proportional to the rate of convection, horizontal temperature variation and the vertical temperature variation. The constants \(\sigma, \rho \), and \( \beta\) are the system parameters proportional to the Prandtl number, Rayleigh number and the physical dimensions of the fluid layer. This system is used to model the lasers, dynamos, thermosyphons, electric circuits, chemical reaction and osmosis. One simple application of these equations is in the modelling of the famous Malkus waterwheel.

We can think of this system as $$\frac{d\textbf{s}}{dt}=\textbf{f}(t,\textbf{s})$$ where \( \textbf{f}(t,\textbf{s})=(\sigma(y-x),x(\rho-z)-y,xy-\beta z) \).

Higher Order ODEs

If we can solve the n-th order differential equaiton for the n-th derivative as follows $$ y^{(n)}=f(t,y, y', y'', y''', \cdots ,y^{(n-1)}), $$ then we can convert the n-th order differential equaiton into a system of n differential equations.

Example 4: Simple Pendulum

Consider the motion of an ideal pendulum with a mass \( m \) attached at the end of the arm of length \( l \). The equation of motion for this pendulum is given As $$ \theta ''+\alpha \theta '+ \beta \sin \theta = \gamma \sin \Omega t,$$ where \( \theta \) is the angle of displacement of the pendulum at tiime \(t \), and \( \alpha, \beta \) and \(\Omega \) are the parameters associated with the system bases on the coefficient of friction, mass of the pendulum, amplitude and the frequency of the forcing term. We will convert the third order differential equaiton in a system of first order equations as follows.

Let $$ x_1 =\theta $$ $$ x_2 = \theta ' $$ Then we have a system of two first order differential equations: $$x_1' = \theta ' =x_2 $$ $$x_2' =\theta '' = -\alpha \theta '- \beta \sin \theta + \gamma \sin (\Omega t)=-\alpha x_2- \beta \sin x_1 + \gamma \sin (\Omega t)$$ The initial conditions are \( \textbf{x}(0) =(x_1(0), x_2(0))^T =(\theta_0, v_0)^T \)

For our numerical computation, we will use the parameters \( \alpha = 0.1, \beta =\gamma = \Omega =1 \). The initial condition is \( (x_1(0), x_2(0)) = (0.5, -2 ) \).

Exercise 3

Solve the following 3rd order nonlinear ODE that arises in a laminar boundary layer problem: $$ff''+f'''=0, \quad f(0)=0, \quad f'(0)=0, \quad f''(0)=0.46960 $$ It is very hard to find an analytical solution to this system, so we will try to solve it numerically.

Disclaimer: THIS PAGE IS NOT A PUBLICATION OF THE UNIVERSITY OF NORTH GEORGIA (UNG) AND UNG HAS NOT EDITED OR EXAMINED THE CONTENT OF THE PAGE. THE AUTHOR OF THE PAGE IS SOLELY RESPONSIBLE FOR THE CONTENT.