author:

• 'Adrian E. Feiguin' title: 'Computational Physics' ...

Ordinary differential equations

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.

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{darray}{rcl} && dx \approx \Delta x &=&x_1-x_0, \\ && dy \approx \Delta y &=&y_1-y_0,\end{darray}$$ 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$.

Exercise 1.1: Newton’s law of cooling

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_s=83^{\circ} C$, $T_0=10^{\circ}$ using different intervals $\Delta t$ and compare with exact (analytical) results.

Let's try plotting the results. We first need to import the required libraries and methods

Next, we create numpy arrays to store the (x,y) values

We have to re write the loop to store the values in the arrays. Remember that numpy arrays start from 0.

We could have saved effort by defining

Alternatively, and in order to re use code in future problems, we could have created a function.

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).

Now, let's study the effects of different time steps on the convergence:

Challenge 1.1

To properly study convergence, one possibility it so look at the result at a given time, for different time steps. Modify the previous program to print the temperature at $t=10$ as a function of $\Delta t$.