. In the Newton's law of cooling example
Euler's method has the form of a first order finite difference algorithm
UCES
Course: Methods and Analysis in UCES
Method: Euler-Trapezoid algorithm
Assessment: second order approximation discretization error
Objectives:
General Information: This section is a continuation of the previous section. Here we focus on stiff problems and a second order accurate discretization scheme.
Contact Person: R. E. White, NCSU, white@math.ncsu.edu
Revision Date: 8-3-94
Copyright ©1994, R. E. White. All rights reserved.
In the previous lecture we fixed the maximum time at T and let the step
length, h, vary. Here we will fix the step length and let T be very large or
. In the Newton's law of cooling example
Euler's method has the form of a first order finite difference algorithm
In the Newton's law of cooling model
In this section we continue to consider a well-stirred liquid so that the temperature is uniform within the liquid. Heat loss to the surrounding region will given by Newton's law of cooling. Here we assume rapid cooling is taking place.
Application to Rapid Cooling. At the end of the
previous section we presented a table for the coffee cup problem where
c = .0162. Both the Euler and Euler-Trapezoid algorithms converged but at the
different predicted rates. Consider
 |
| Figure: Euler Algorithm: Blowup with c = .405 |
The model is from Newton's law of cooling. Let u(t) be the temperature of the liquid.
The constant c determines how fast heat is leaving the container. For an insulated cup c will be relatively small. If the cup is made of metal, c will be larger. Further, if the cup is metal and there is a large breeze, then c will be even larger.
The Euler-Trapezoid algorithm, which can compute the
solution of stiff differential equations with second order accuracy, is based
on the trapezoid rule for numerical integration. Note,
The discrete form of this is
Observe the right side where the unknown uk+1 also appears. This is an additional complication relative to Euler's algorithm. In a previous section we studied the bisection algorithm which was used to find roots. There are two other root finding methods which we will study in the next two sections. For the present we will use a Maple procedure called fsolve to find the unknown.
There are two important differences between the Euler algorithm and the Euler-Trapezoid algorithm. First, the Euler algorithm has an error bound which is proportional to h; where as, the Euler-Trapezoid algorithm has an error bound which is proportional to h2 . Second, in the case that fu is negative the error bound for the Euler-Trapezoid algorithm is valid for all k and not just those between 1 and K with h = T/K. This is particularly important when fu is a large negative number.
Euler-Trapezoid Algorithm.
choose h
u0 = u(0)
for k = 1,maxk
oldf = f(k*h,u(k) )
solve for newu = u(k) + (h/2)*(oldf + f((k+1)*h,newu))
u(k+1) = newu
endloop.
This is a modified version of the Maple code for the Euler algorithm to the
Euler-Trapezoid algorithm. We have used the Maple procedure fsolve to solve
the nonlinear problem associated with each time step. The first graph is for
the same nonstiff data that was used in the Euler implementation in the
previous section. We have also used a variety of mesh sizes. Note the
Euler-Trapezoid algorithm has converged much more rapidly than the Euler
algorithm.
Maple Code for the Euler-Trapezoid Algorithm.
Next we try some stiff data that was reported earlier where c = .405
and with variable mesh h = 5 (K = 10), 2.5 (K = 20), and 1.25 (K = 40). Later
we shall use an alternate method for finding the fixed point. .
The next computations are done using Matlab, and we illustrate the more rapid
convergence of the Euler-Trapezoid algorithm for the differential equation
from Newton's law of cooling where
Solve for
The numerical errors for different numbers of steps, kk (h = dt = T/kk), are
computed at time equal to 1.0, and the reader should compare the following
figure of the Euler-Trapezoid error with figure in the previous section for the
Euler error.
Matlab Code for Newton's Law of Cooling.
We now give the discretization error analysis for the special case of Newton's
law of cooling where
Now use the mean value theorem
Solve for u(kh) and insert it into the above.
At this point we again use the backward substitution
Assume
Euler-Trapezoid Error Theorem. Consider the
Euler-Trapezoid algorithm as it is applied to
If the first and second derivatives are bounded, and M and a are defined as in
(1), then (2) must hold for all k.
2.5.5. Implementation.
Input Data:
> t:='t':
> u:='u':
> f:=(t,u)->.061*(70 - u);
> t0:=0:
> u0:=200:
> T:=50:
> K:=10:
Procedure EULERT is Defined:
> EULERT:=proc(f,t0,u0,T,K)
> t(0):= t0:
> h:=T/K:
> u(0):= u0:
> for k from 0 to K do
> oldf:=f(t(k),u(k)):
> t(k+1):= t(k) + h:
> fsolve(x=u(k) + (h/2)*(oldf + f(t(k+1),x)),x):
> u(k+1):=x:
> od:
> end:
Output Data:
> with(plots):
> EULERT(f,0,200,50,4):
> L:=[seq([t(n),u(n)],n=0..4)]:
> EULERT(f,0,200,50,8):
> LL:=[seq([t(n),u(n)],n=0..8)]:
> EULERT(f,0,200,50,16):
> LLL:=[seq([t(n),u(n)],n=0..16)]:
> plot({L,LL,LLL},x=0..50);
Figure: Euler-Trapezoid Solutions
Figure: Solution of Stiff Equation
newu
= u(k) + (h/2)*(oldf + f((k+1)*h,newu))
= u(k) + (h/2)*(c*(usur - u(k)) + c*(usur - newu)).
clear;
i=0;
vec = zeros(12,1);
vect = 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;
ut = u0;
for k = 1:kk
uk = uk +dt*c*(usur -uk);
uex = usur + (u0 -usur)*exp(-c*k*dt);
ut = (ut +dt*c*usur - dt*c*ut/2)/(1+c*dt/2);
error = abs(uk - uex);
errort = abs(ut - uex);
end
i = i+1;
vect(i) = errort;
vec(i) = error;
vecdt(i) = 1/dt;
end
plot(vecdt,vect)
Figure: Euler-Trapezoid Error versus Steps
.
(1)
and use the fact that
a < 1 to conclude
(2)
.