UCES

"Initial Value Problems: Euler's Method"

Course: Methods and Analysis in UCES

Chapter:   2.0   Computational Calculus
Section:   2.4   Initial Value Problems: Euler's Method
Introduction:   initial value problem

Applied Area:   heat transfer

Model:   Newton's law of cooling

Method:   Euler's algorithm

Implementation:   Maple, Matlab

Assessment:   discrete error, accumulation errors

Homework:

Prerequisites:

1. Math:   Calculus
2. Computing:   Maple and loops
3. Science:   Physics

Objectives:

1. Math:   Convergence analysis and extended mean value theorem Computing:   Use of loops Science:   First continuum model of heat transfer

General Information: This is the first of a series of continuum models for heat transfer. We will illustrate model evolution and the approximation of continuum models by discrete models.

Contact Person:   R. E. White, NCSU, white@math.ncsu.edu

Revision Date:   8-2-94

Copyright ©1994, R. E. White. All rights reserved.

Return to Table of Contents

2.4   Initial Value Problems: Euler's Method

Initial value problems have the form

ut = f(t,u) and u(0) given.

(1)

The simplest cases can be solved by separation of variable, but in general they do not have to have closed form solutions; for example, let f(t,u) = t² + u². Therefore, one is forced to consider various approximation methods. In the next two lectures we present the simplest numerical methods. Both methods are based upon approximating the initial value problems by a sequence of calculations which divide to time interval [0,T] into K parts. The mesh size is h = T/K. We shall see that the error made by this type of approximation is proportional to h for the Euler's method, and to h² for the second method, Euler-Trapezoid method.

2.4.2.   Applied Area.

Consider a well stirred liquid such as a cup of coffee. We will assume that the temperature is uniform with respect to space, but the temperature may be changing with respect to time. We wish to predict the temperature as a function of time given some initial observations about the temperature. For example, a cup of coffee will cool much faster in a metal cup than in an insulated cup. Or, the coffee will cool much faster if the surrounding temperature is very cold.

2.4.3.   Model

Newton's empirical law of cooling states that the rate of change of the temperature is proportional to the difference in the surrounding temperature and the temperature of the liquid. This is a continuum version of the one described in the first section on floating point numbers. In fact, the Euler algorithm for the following differential equation is exactly the discrete version!

u_t = c(u_sur - u) where

u is the temperature of the liquid,
u(0) is observed,
u_sur is the surrounding temperature and
c is the coefficient which reflects the insulation of the container.

The closed form solution of this differential equation can be found by the separation of variables method and is

u(t) = u_sur + (u(0) - u_sur )e^-ct .

If the c is not given, then it can be found from a second observation u(t₁ ) = u₁ . In this case one can approximate c by a finite difference approximation which we will later illustrate. If u_sur is a function of t, one can still find a closed form solution provided it is not too complicated to be easily integrated.

2.4.4.   Method.

The first method involves the approximation of u_t by the finite difference

(u^k+1 - u^k )/h

where h = T/K, u^k is an approximation of u(kh) and f is evaluated at (kh, u^k ). If T is not finite, then h will be fixed and k may range over all of the positive integers. The differential equation (1) can be replaced by either

(u^k+1 - u^k )/h = f((k+1)h, u^k+1 )

(uk+1 - uk)/h = f(kh,uk).

(2)

The second choice is the simplest because it does not require the solution of a possibly nonlinear problem at each time step. The scheme given by (2) is called Euler's method.

Let u(k) be an array of dimension K+1 with u(0) equal to the initial temperature.

u(k) will be the approximate values of the temperature at time equal to kh , that is, u(k) = u^k . Euler's algorithm can be viewed as a generalization of the first order finite difference algorithm where u + hf(t,u) = au + b, or solving for f, f(t,u) = ((a-1)u + b)/h.

Euler's Algorithm.

        h = T/K
	u(0) = u0
	for k = 0, K-1
		u(k+1) = u(k) + h*f(k*h,u(k))
	endloop.

2.4.5.   Implementation.

Use Newton's law of cooling to predict the temperature of a cup of coffee. The coffee is kept well stirred by walking around the classroom. Also, we need to know when its temperature reaches 110, after that the professor gets irritable. Today the initial temperate was 200 and five minutes later it was observed to be 190. The second measurement varies with the type of coffee container. The continuum model is

u_t = c(u_sur - u) where

u(0) = 200, u(5) = 190 and u_sur = 70.

It is possible to solve this model in closed form, and the solution is

u(t) = 70 + (200 - 70)e^-ct .

The c is determined from the second observation and is c = (-1/5)ln(12/13) = .01600854

The discrete model is formed by using Euler's algorithm and is

u^k+1 = u^k + hc(u_sur - u^k) where

u0 = 200	usur = 70
T = 50	h = 50/K.

The c must be determined from the second observation. Here the increase of time from 0 to 5 is not too large, and so we may approximate the u_t (0) by (190 - 200)/5 and now put this into the differential equation

(190 - 200)/5 = c(70 - 200) or c(71 - 190) or c((70 + 71)/2 - (200 + 190)/2).

We will choose the last possible right side so that c = 2/123.5 = .01619433 which is not too far from the c in the continuum model. For K =10 the time step is h = 50/10, and hc = .08097. At the end of the class the coffee's computed temperature was 127, yeah.

The computations below use a larger c = .061 and indicate the effect of the mesh size on the numerical solution. For this example the numerical solution appears to converge as K increases from 4 to 8 and to 16. The graphs of the numerical solution are in increasing order.

Maple Code for a Well Stirred Liquid.

  The Procedure is Defined:
	> EULER:=proc(f,t0,u0,T,K)
	> t(0):= t0:
	> h:=T/K:
	> u(0):= u0:
	> for k from 0 to K do
	> 	u(k+1):= u(k) + h*f(t(k),u(k)):
	> 	t(k+1):= t(k) + h:
	> od:
	> end:

  Input Data:
	> t:='t':
	> u:='u':
	> f:=(t,u)->.061*(70 - u);
	> t0:=0:
	> u0:=200:
	> T:=50:
	> K:=10:

The Procedure is Called for Different Mesh h = T/K:
	> with(plots):

	Use K = 4.
	> EULER(f,0,200,50,4):
	> L:=[seq([t(n),u(n)],n=0..4)]:

	Use K = 8.
	> EULER(f,0,200,50,8);
	> LL:=[seq([t(n),u(n)],n=0..8)]:

	Use K = 16.
	> EULER(f,0,200,50,16);
	> LLL:=[seq([t(n),u(n)],n=0..16)]:

  The Output is Graphed for the Three Mesh Choices:
	> plot({L,LL,LLL},x=0..50);

Figure: Convergence of Numerical Solution

The next numerical experiments are done using Matlab. They illustrate numerical error as a function of the number of steps. This is done with the equation for Newton's law of cooling with usur = 70 and c = 1.0; the error is computed for time equal to 1.0.

Matlab Code for Euler's Algorithm.

	clear;
	i=0;
	vec = zeros(12,1);
	vecdt= zeros(12,1);
	for kk=4:4:48
		T  =1.0;
		dt =T/kk;
		u0 = 200.;
		c  = 1.0;
		usur  = 70.;
		uk  = u0;
		for k = 1:kk
			uk = uk +dt*c*(usur -uk);
			uex = usur + (u0 -usur)*exp(-c*k*dt);
			error = abs(uk - uex);
		end
		i= i+1;
		vec(i) = error;
		vecdt(i) = 1/dt;
	end
	plot(vecdt,vec)

Figure: Error versus Steps (h = T/steps)

2.4.6.   Assessment.

If the liquid is not well stirred, then the temperature will vary within the container. This results in diffusion of heat from a warm region to a cool region. The above model does not take this into consideration. Also, in our analysis the surrounding temperature was held constant. Variable surrounding temperature can easily be accounted for in our present model.

Clearly this is a very inexpensive algorithm to execute. However, there are sizable errors. There are two types of errors: the discretization and accumulation errors.

Discretization error = E_d^k = u^k - u(kh) where

u^k is from Euler's algorithm (2) with no roundoff error and

u(kh) is from the exact continuum solution (1).

Accumulation error = E_r^k = U^k - u^k where

U^k is from Euler's algorithm, but with roundoff errors.

2.4.7.   Homework.

1. Fill in the steps leading to (6).

2. Use (5) and the arguments similar to those leading to the discretization error, to derive the accumulation error estimate in the theorem.

3. Write a computer code for Euler's algorithm, use the above example with

   f(t,u) = t² + u².

   Experiment with different initial conditions, u⁰, and different step lengths, h.

4. Consider the above coffee cup problem, but now the classroom temperature is observed to increase in a linear way from 70 to 80 at the end of class 50 minutes later. Find both the continuum and discrete solutions, and do the computations similar to those given in the above table.