- A is the area , say A = 1,
= .022 is the emissivity,
= 5.68 10-8 is the Stefan-Boltzmann
constant and
- usur = 273 is the surrounding temperature.
UCES
Course: Methods and Analysis in UCES
Prerequisites:
Objectives:
General Information: This and the next section present alternatives to the bisection algorithm for finding roots and fixed points. The Picard and Newton methods generalize in many ways to a variety of similar problems.
Contact Person: R. E. White, NCSU, white@math.ncsu.edu
Revision Date: 8-2-94
Copyright ©1994, R. E. White. All rights reserved.
The Euler-Trapezoid algorithm requires finding at each time step a fixed point of a function x = g(x), or equivalently, a root of f(x) = x - g(x) = 0. This is a common problem that arises in computations, and a more general problem is to find N unknowns when N equations are given. The bisection algorithm does not generalize very well to these more complicated problems. In this and the next section we will present two algorithms, Picard and Newton, which do generalize.
In the previous section we considered the rapid cooling of an object which had uniform temperature with respect to the space variable. Heat loss by transfer from the object to the surrounding region may be governed by equations that are different from Newton's law.
Application to Radiative Cooling. Suppose a thin plate is glowing hot so that the main heat loss is via radiation. Then Newton's law of cooling may not be an accurate model. A more accurate model is the Stefan radiation law
= .022 is the emissivity,
= 5.68 10-8 is the Stefan-Boltzmann
constant and
The derivative of f is -4cu3 and is large and negative for temperature near the initial temperature, f'(973) = -4.6043; but it is smaller for temperatures near the surrounding temperature, f'(273) = -0.1017. Consequently, it is initially a stiff differential equation, but as the temperature approaches the surrounding temperature, it is less stiff.
We will need to use an algorithm which is suitable for stiff differential equations. The Euler-Trapezoid algorithm is suitable provided we can solve the associated fixed point problem at each time step.
 |
(1) |
For small enough h this can be solved by the Picard algorithm (more on this in the next section)
| um+1 = g(um ) | (2) |
where the m indicates an inner iteration and not the time step. The initial guess for this iteration can be taken from one step of the Euler algorithm.
Example. Consider the first time step for
This can be solved using the quadratic formula, but for small enough h one can
use (2) several times. Let h = .1 and let the first guess be u0 = 1
(m = 0). Then the calculations from (2) will be: 1.100500, 1.111055,
1.112222, 1.112351. If we are "satisfied" with the last calculation, we let it
be the value of the next time set, uk where
Consider the problem of finding the fixed point of g(x) = x.
The Picard algorithm has the form of successive approximation as in (2), but for more general g(x). This algorithm is also called the contractive algorithm for reasons that we will soon see. In the algorithm below we continue to iterate (2) until there is little difference, eps, in two successive calculations. The eps is positive and "suitably" small.
choose oldx in [a,b] and eps > 0
for m = 1, maxit
newx = g(oldx)
if |xnew - xold| < eps go to end
oldx = newx
endloop
end.
Example. Find the square root of 2. This could also be written either as 0 = 2 - x2 for the root, or as x = x + 2 - x2 = g(x) for the fixed point of g(x). Try an initial approximation of the fixed point, say, x0 = 1. Then the subsequent iterations are
So, the iteration does not converge. Try another initial x0 = 1.5 and it still does not converge
Note, this last sequence of numbers is diverging from the solution. A good way to analyze this is to again use the mean value theorem
) = g'(c)(x - )
.
Here g'(c) = 1 - 2c. So, regardless of how close x is to
, g'(c) will approach
we have
)| >
|x - |!
Definition. g:[a,b]
R
is called contractive on [a,b] if and only if for all x and y in [a,b]
and positive
| g(x) - g(y) |  r | x - y |. |
(3) |
Example. Consider a simple implicit discretization of
One can verify that the first 6 iterates of Picard's algorithm with h = 1 are
The algorithm has converged to within 10-4, and we stop and set
In this problem for all values of
1
Euler-Trapezoid Algorithm with Picard Solver.
choose h and eps u0 = u(0) for k = 1,maxk oldf = f(k*h,u(k) ) oldu = u(k) + h*oldf for m = 1, maxit newu = u(k) + (h/2)*(oldf + f((k+1)*h,oldu)) if |oldu - newu| < eps goto end oldu = newu endloop end u(k+1) = newu endloop.
Example. Let us return to the radiative cooling problem where we must
solve a stiff differential equation. If we use the above algorithm with a
Picard solver, then we will have to find a fixed point of
max |g'(c)| |x - y|.
So, we must require r = max |g'(c)| < 1. In our case,
This is a modified version of the Maple code for the Euler algorithm to the Euler-Trapezoid algorithm. Note the inner loop which solves the fixed point problem by Picard iteration. The first graph is for the same nonstiff data that was used in the Euler implementation in the previous lecture. 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 with Picard Solver.
Input Data:
> t:='t':
> u:='u':
> f:=(t,u)->.6941(1 - (u/273)^4)):
> t0:=0:
> u0:=973:
> T:=100:
> K:=10:
Procedure EULERTP is Defined:
> EULERT:=proc(f,t0,u0,T,K)
> t(0):= t0:
> h:=T/K:
> u(0):= u0:
> eps:=.0001:
> maxm:= 50:
> for k from 0 to K do
> oldf:=f(t(k),u(k)):
> oldu:=u(k) + h*oldf:
> t(k+1):= t(k) + h:
> test:=1:
> for m from 1 to maxm while(test > eps) do
> newu:=u(k) + (h/2)*(oldf + f(t(k+1),oldu)):
> test:= abs(oldu - newu):
> oldu:=newu:
> od:
> u(k+1) :=newu:
> od:
> end:
Output of Data:
> EULERT(f,0,973,100,10):
> L:=[seq([t(n),u(n)],n=0..10)]:
> EULERT(f,0,973,100,20):
> LL:=[seq([t(n),u(n)],n=0..20)]:
> EULERT(f,0,973,100,40):
> LLL:=[seq([t(n),u(n)],n=0..40)]:
> plot({L,LL,LLL},x=0..50);
 |
| Figure: Temperature versus Time: K = 10, 20, 40 |
Picard's algorithm converged even for the largest values of h = 10. The theory requires h to be less than one half for u near 973. As the temperature decreases, the step size h does not have to be as small.
In the radiative cooling model we have also ignored the good possibility that there will be differences in temperature according to the location in space. In such cases there will be diffusion of heat, and one must model this mode of heat transfer.
We indicated that the Picard algorithm will converge if the mapping g(x) is contractive. The following theorem makes this more precise. The error is order one convergent, provided g is contractive and its derivative is bounded. This means the new error is bounded by the old error
r | xm - x |.
Picard Convergence (Contraction Mapping) Theorem.
Let g:[a,b]
[a,b] and assume that x is a fixed point
of g and x is in [a,b]. If g is contractive on [a,b], then the Picard
algorithm converges to the fixed point. Moreover, the fixed point is unique.
Proof. Let xm+1 = g(xm ) and x = g(x). Repeatedly use the contraction property (3).
| | xm+1 - x | | = | g(xm ) - g(x) | |
| r | xm - x | |
|
| = r | g(xm-1 ) - g(x) | | |
| r2 | xm-2
- x | |
|
| . . . |
|
| rm+1 | x0
- x |. |
r < 1, rm+1 must go to
zero as m increases.
If there is a second fixed point, y, then
r | y - x
| r