Objectives:
- Math: Newton's method in more than one variable
- Computing: Solution of nonlinear problem
- Science: Nonlinear thermal properties
General Information: This is an extension of the one variable Newton method for finding roots of a single equation to finding solutions for n variables with n equations. This very important method is applied to steady state 1D space diffusion with nonlinear cooling due to radiation.
Contact Person: R. E. White, NCSU, white@math.ncsu.edu
Revision Date: 8-4-94
Copyright ©1994, R. E. White. All rights reserved.
3.7 Nonlinear Boundary Value Problems and Radiative Heat Transfer
3.7.1. Introduction.
In the analysis of most of the heat transfer problems we assumed the temperature varied over a small range so that the thermal properties could be approximated by constants. This always resulted in a linear algebraic problem which could be solved by a variety of methods. Two possible difficulties are (i) nonlinear thermal properties, or (ii) larger problems which are a result of diffusion in three directions. In this section we consider the nonlinear problems.
The thermal properties of specific heat and thermal conductivity can be nonlinear. The exact nature of the nonlinearity will depend on the material and the range of the temperature variation. Usually, data is collected that reflects these properties, and a least squares curve fit is done for a suitable approximating function. Other nonlinear terms can evolve from the heat source or sink terms in either the boundary conditions, or the source term on the right side of the heat equation. We consider one such case.
Consider a cooling fin or plate which is glowing hot, say at 900 degrees Kelvin. Here heat is being lost by radiation to the surrounding region. In this case the heat lost is not proportional, as in Newton's law of cooling, to the difference in the surrounding temperature, usur , and the temperature of the glowing mass, u. Observations indicate that the heat loss through a surface area, A, in a time interval, dt, is equal to
- dt A
(usur4 - u4 )
is the emissivity of the surface and
is the
Stefan-Boltzman constant. One can couple this with the Fourier heat law to
form various differential equations or boundary conditions which are
nonlinear.
Radiative Cooling in a Thin Wire. Consider a
thin wire of length L and radius r. Let the ends of the wire have a fixed
temperature of 900 and let the surrounding region be
r dx.
Therefore, the heat leaving the lateral surface in dt time is
- dt (2r dx) ((usur4 -
u4 )) .
Assume steady state heat diffusion in one direction and apply the Fourier heat law to get
- 0
dt (2r dx)
(
(usur4 - u4 )) +
dt K(r2 )ux(x + dx/2) - dt
K(r2 )ux(x - dx/2) .
Now, divide by dt dx r and let dx go to zero to get
The continuum model for the heat transfer is
The thermal conductivity will also be temperature dependent, but we will
assume for simplicity that K is a constant which has been incorporated into
c.
Consider the nonlinear differential equation
This has n unknowns, ui , and n equations of the form
Find u = [u1 , ..., un ]T so that
Nonlinear problems can have multiple solutions. For example, consider the
intersection the unit circle
In order to derive Newton's method for n equations and n unknowns, it is
instructive to review the one equation case. The idea behind Newton's
algorithm is to approximate the function f(x) at a given point by a straight
line. Then find the root of the equation associated with this straight line.
One continues to repeat this until little change in the approximations of the
root is observed. We make use of the approximation
The equation for the straight line at the mth iteration is
Define xm+1 so the y = 0 where the straight line intersects the x
axis
Solve for xm+1 to obtain Newton's algorithm
The derivation of Newton's method for more than one equation and one unknown
requires an analog of the approximation in (3). Consider Fi (u) as
a function of n variables uj. If just the j component of u changes,
then (3) will hold with x replaced by uj and f replaced by
Fi. If all of the components change, then the net change in
Fi can be approximated by sum of the (partial derivatives of
Fi with respect to uj ) times the (change in
uj ):
For n = 2 this is depicted in the following figure where
This can be put into matrix form to obtain the desired analog of (3)
Newton's method is obtained by letting
The vector approximation in (4) is replaced by an equality to get
This vector equation can be solved for um+1 , and we have the
n variable Newton method
In practice the inverse of the Jacobian in not used, but we must find the
solution,
Consequently, Newton's method consists of solving a sequence of linear
problems. One usually stops when either F is "small", or
Example. Let n = 2,
and it is nonsingular if both variables are nonzero.
If the initial guess is near the solution in a quadrant, then Newton's method
will converge to that solution.
Example. Consider the nonlinear differential equation for the radiative
heat transfer problem in (2) where
The Jacobian matrix is easily computed and must be tridiagonal because each
Fi only depends on
The following is a Maple code which uses Newton's method to solve the 1D
diffusion problem with heat loss due to radiation. We have used the Maple
procedure linsolve to solve each linear subproblem. One could use
an iterative method, and for a larger problem where there is diffusion of heat
in more than one direction this might be the best way. We have experimented
with several emissitivities, and the graphs indicate the higher the
emissitivity the more the cooling.
Maple Code for Nonlinear Cooling with Diffusion.
The next calculations were done using the Matlab version of the above
code, and the very rapid convergence of Newton's method is illustrated.
Matlab Code for Nonlinear Cooling with Diffusion.
Nonlinear problems are very common, but they are often linearized. This is
done because it is easier to solve one linear problem than a nonlinear
problem. In the case of Newton's method we must solve a sequence on linear
problems. If one tries to use analytic methods for finding closed form
solutions of nonlinear boundary value problems, it may be impossible!
Newton's method has, under some assumption on F(u), the two very important
properties of local convergence and quadratic convergence.
Local convergence means that if the initial
guess is suitably close to the solution, then the Newton iterates will
converge to the solution. Often one can use physical insight to the problem
to make a good initial guess. Also, such initial guesses can significantly
reduce the number of Newton iterations needed to achieve convergence.
Quadratic convergence means the error at the
next step will be proportional to the square of error of the present step.
Thus, after the error is less than 1, the convergence is very rapid. These
two properties have contributed to the wide use and many variations of the
Newton's method for solving nonlinear algebraic systems. In order for the
above properties to hold the components functions should be sufficiently
differentiable, and the Jacobian matrix should be inevitable and well
conditioned.
Use Newton's method coupled with a linear solver that uses SOR to solve
this nonlinear problem.
(usur4 - u4 ) +
(Krux )x.
Figure: Wire Segment
-(Kux )x = c(usur4 -
u4 ) where c = 2
/r and(1)
u(0) = 900 = u(L) .
-ui-1 + 2ui - ui+1 =
h2f(ui ) for i = 1, ..., n.
(2)
df
f'(x) dx.(3)
Figure: Newton's Algorithm
Figure: Change in Fi
where
F
F'(u)
u(4)
F =
[F1 , ...,
Fn]T,u = [u1 , ..., un]T
u = um+1 - um and
F = 0 - F(um ) .
u, of
u.
u is "small."
choose initial u
for m = 1,maxit
compute F(u) and F'(u)
solve F'(u)
u =-F(u)
u = u + u
test for convergence
endloop.
Input Data:
> with(linalg):
> with(plots):
> f:=(u)->.000005*(300^4 - u^4);
> fp:=(u)->.000005*(-4*u^3);
> n:=19:
> h:=1/(n+1):
> u0:=900:
> eps:=.1:
Procedure NEWTON is Defined:
> NEWTON:=proc(n,h,eps,f,fp,u0):
> FP:=matrix(n,n,0):
> F:=vector(n,0):
> u:=vector(n,u0):
> for m from 1 to 10 do
> for i from 1 to n do
> if i=1 then
> F[i]:=f(u[i])*h*h - 2*u[i] + u[i+1] + u0:
> FP[i,i]:=fp(u[i])*h*h -2:
> FP[i,i+1]:=1:
> fi:
> if i>1 and i<n then
> F[i]:=f(u[i])*h*h + u[i-1] - 2*u[i] + u[i+1]:
> FP[i,i]:=fp(u[i])*h*h - 2:
> FP[i,i-1]:=1:
> FP[i,i+1]:=1:
> fi:
> if i=n then
> F[i]:=f(u[i])*h*h - 2*u[i] + u[i-1] + u0:
> FP[i,i]:=fp(u[i])*h*h -2:
> FP[i,i-1]:=1:
> fi:
> od:
> du:=linsolve(FP,F):
> u:=evalm(u - du):
> error:=norm(F):
> if error<eps then break fi:
> od:
> end:
NEWTON is Called and Output of Data:
> f:=(u)->.000005*(300^4 - u^4):
> fp:=(u)->.000005*(-4*u^3):
> NEWTON(n,h,eps,f,fp,u0):
> m;
8
> L:=[[0,u0],seq([i*h,u[i]],i=1..n),[1,u0]]:
> f:=(u)->.0000001*(300^4 - u^4):
> fp:=(u)->.0000001*(-4*u^3):
> NEWTON(n,h,eps,f,fp,u0):
> m;
6
> LL:=[[0,u0],seq([i*h,u[i]],i=1..n),[1,u0]]:
> f:=(u)->.00000005*(300^4 - u^4):
> fp:=(u)->.00000005*(-4*u^3):
> NEWTON(n,h,eps,f,fp,u0):
> m;
5
> LLL:=[[0,u0],seq([i*h,u[i]],i=1..n),[1,u0]]:
> plot({L,LL,LLL},x=0..1);
Figure: Temperatures for Different
Emissitivities
uo = 900.;
n = 19;
h = 1./(n+1);
FP = zeros(n);
F = zeros(n,1);
u = ones(n,1)*uo;
for m =1:20
for i = 1:n
if i==1
F(i) = f(u(i))*h*h + u(i+1) - 2*u(i) + uo;
FP(i,i) = fp(u(i))*h*h - 2;
FP(i,i+1) = 1;
elseif i<n
F(i) = f(u(i))*h*h + u(i+1) - 2*u(i) + u(i-1);
FP(i,i) = fp(u(i))*h*h - 2;
FP(i,i-1) = 1;
FP(i,i+1) = 1;
else
F(i) = f(u(i))*h*h - 2*u(i) + u(i-1) + uo;
FP(i,i) = fp(u(i))*h*h - 2;
FP(i,i-1) = 1;
end
end
du = FP\F;
u = u - du;
error = norm(F);
if error<.0001
break;
end
m
error
end
Output Data:
m = 1
error = 706.1416
m = 2
error = 197.4837
m = 3
error = 49.2847
m = 4
error = 8.2123
m = 5
error = .3967
m = 6
error = .0011
m = 7
error = 7.3703e-09