author:
Let’s consider a simple 1st order equation: $$\frac{dy}{dx}=f(x,y)$$
To solve this equation with a computer we need to discretize the differences: we have to convert the differential equation into a “finite differences” equation. The simplest solution is Euler’s method.
Supouse that at a point $x_0$, the function $f$ has a value $y_0$. We want to find the approximate value of $y$ in a point $x_1$ close to $x_0$, $x_1=x_0+\Delta x$, with $\Delta x$ small. We assume that $f$, the rate of change of $y$, is constant in this interval $\Delta x$. Therefore we find: $$\begin{eqnarray} && dx \approx \Delta x &=&x_1-x_0, \\ && dy \approx \Delta y &=&y_1-y_0,\end{eqnarray}$$ with $y_1=y(x_1)=y(x_0+\Delta x)$. Then we re-write the differential equation in terms of discrete differences as: $$\frac{\Delta y}{\Delta x}=f(x,y)$$ or $$\Delta y = f(x,y)\Delta x$$ and approximate the value of $y_1$ as $$y_1=y_0+f(x_0,y_0)(x_1-x_0)$$ We can generalize this formula to find the value of $y$ at $x_2=x_1+\Delta x$ as $$y_{2}=y_1+f(x_1,y_1)\Delta x,$$ or in the general case: $$y_{n+1}=y_n+f(x_n,y_n)\Delta x$$
This is a good approximation as long as $\Delta x$ is “small”. What is small? Depends on the problem, but it is basically defined by the “rate of change”, or “smoothness” of $f$. $f(x)$ has to behave smoothly and without rapid variations in the interval $\Delta x$.
Notice that Euler’s method is equivalent to a 1st order Taylor expansion about the point $x_0$. The “local error” calculating $x_1$ is then $O(\Delta x^2)$. If we use the method $N$ times to calculate $N$ consecutive points, the propagated “global” error will be $NO(\Delta x^2)\approx O(\Delta x)$. This error decreases linearly with decreasing step, so we need to halve the step size to reduce the error in half. The numerical work for each step consists of a single evaluation of $f$.
If the temperature difference between an object and its surroundings is small, the rate of change of the temperature of the object is proportional to the temperature difference: $$\frac{dT}{dt}=-r(T-T_s),$$ where $T$ is the temperature of the body, $T_s$ is the temperature of the environment, and $r$ is a “cooling constant” that depends on the heat transfer mechanism, the contact area with the environment and the thermal properties of the body. The minus sign appears because if $T>T_s$, the temperature must decrease.
Write a program to calculate the temperature of a body at a time $t$, given the cooling constant $r$ and the temperature of the body at time $t=0$. Plot the results for $r=0.1\frac{1}{min}$; $T_0=83^{\circ} C$ using different intervals $\Delta t$ and compare with exact (analytical) results.
In [2]:
T0 = 10. # initial temperature
Ts = 83. # temp. of the environment
r = 0.1 # cooling rate
dt = 0.05 # time step
tmax = 60. # maximum time
nsteps = int(tmax/dt) # number of steps
In [3]:
T = T0
for i in range(1,nsteps+1):
new_T = T - r*(T-Ts)*dt
T = new_T
print i,i*dt, T
# we can also do t = t - r*(t-ts)*dt
Let's try plotting the results. We first need to import the required libraries and methods
In [6]:
%matplotlib inline
import numpy as np
from matplotlib import pyplot
Next, we create numpy arrays to store the (x,y) values
In [7]:
my_time = np.zeros(nsteps)
my_temp = np.zeros(nsteps)
We have to re write the loop to store the values in the arrays. Remember that numpy arrays start from 0.
In [8]:
T = T0
my_temp[0] = T0
for i in range(1,nsteps):
T = T - r*(T-Ts)*dt
my_time[i] = i*dt
my_temp[i] = T
In [79]:
pyplot.plot(my_time, my_temp, color='#003366', ls='-', lw=3)
pyplot.xlabel('time')
pyplot.ylabel('temperature');
We could have saved effort by defining
In [80]:
my_time = np.linspace(0.,tmax,nsteps)
In [81]:
pyplot.plot(my_time, my_temp, color='#003366', ls='-', lw=3)
pyplot.xlabel('time')
pyplot.ylabel('temperature');
Alternatively, and in order to re use code in future problems, we could have created a function.
In [9]:
def euler(y, f, dx):
"""Computes y_new = y + f*dx
Parameters
----------
y : float
old value of y_n at x_n
f : float
first derivative f(x,y) evaluated at (x_n,y_n)
dx : float
x step
"""
return y + f*dx
In [11]:
T = T0
for i in range(1,nsteps):
T = euler(T, -r*(T-Ts), dt)
my_temp[i] = T
Actually, for this particularly simple case, calling a function may introduce unecessary overhead, but it is a an example that we will find useful for future applications. For a simple function like this we could have used a "lambda" function (more about lambda functions here).
In [12]:
euler = lambda y, f, dx: y + f*dx
Now, let's study the effects of different time steps on the convergence:
In [21]:
dt = 1.
#my_color = ['#003366','#663300','#660033','#330066']
my_color = ['red', 'green', 'blue', 'black']
for j in range(0,4):
nsteps = int(tmax/dt) #the arrays will have different size for different time steps
my_time = np.linspace(dt,tmax,nsteps)
my_temp = np.zeros(nsteps)
T = T0
for i in range(1,nsteps):
T = euler(T, -r*(T-Ts), dt)
my_temp[i] = T
pyplot.plot(my_time, my_temp, color=my_color[j], ls='-', lw=3)
dt = dt/2.
pyplot.xlabel('time');
pyplot.ylabel('temperature');
pyplot.xlim(8,10)
pyplot.ylim(48,58);