Parallel Spectral Numerical Methods/Nonlinear Ordinary Differential Equations and Iteration

Nonlinear Ordinary Differential Equations and IterationEdit

The implicit explicit method avoids the direct solution of nonlinear problems. This can be advantageous for some problems, but can also lead to severe time step restrictions in others. Furthermore, the resulting numerical schemes can sometimes have undesirable qualitative properties. For this reason, we need to describe methods that allow us to solve the nonlinear equations generated in fully-implicit numerical schemes.

We consider an ordinary differential equation   for  , and for which   is not necessarily a linear function of  . We want to use an implicit numerical method to obtain an approximate solution of this problem – for example backward Euler’s method. If we want to demonstrate the convergence of the numerical scheme, we need to demonstrate convergence of functional iteration which we use to find the solution for the nonlinear equation term in using backward Euler’s method.

The results that follow are primarily taken from Iserles[1], although this material is also often found in calculus texts such as Lax, Burstein and Lax[2], and Hughes et al.[3]. We will let   denote the time at time step  ,   denote the approximate solution at time step   and   denote the time step. We will assume   is Lipschitz continuous, a condition that is weaker than differentiable but stronger than continuous, which we will give a precise definition of. There are two classical iteration methods:

  • fixed-point iteration
  • Newton’s (Newton-Raphson) method.

We will prove convergence of these two methods (a proof of the convergence of the modified Newton-Raphson method is in Iserles[4]). We will analyze the specific problem   with initial data   and  .

Exact Solution to an Example Nonlinear Ordinary Differential EquationEdit

We consider   with initial data   and  . Whenever the solution   exists, it will be positive all the time, because the initial value is positive and   is positive.

To integrate this equation explicitly, we use separation of variables to find that   which implies   where   is the constant of integration. Using our initial data we get  , so   is our exact solution for this problem. We will use this exact solution to compare the numerical solutions obtained by the different iterative methods. Notice that this exact solution becomes infinite as  .

Definitions Required to Prove ConvergenceEdit

DefinitionA function   is Lipschitz if   for all   and   in the domain  .

There are two specific definitions of the Lipschitz condition.

Definition The function   is called locally Lipschitz if, for each  , there exists an   such that   is Lipschitz on the open ball of center   and radius  .

DefinitionIf   is Lipschitz on all of the space   (i.e. The open ball is   in above definition), then   is globally Lipschitz.

Note the fundamental difference between the local and global versions of the Lipschitz-condition. Whereas in the local version the Lipschitz “constant”(  and the open ball depend on each point   , in the global version the “constant”(  is fixed and the open ball is  . In particular, a globally Lipschitz function is locally Lipschitz continuous, but the converse is not true.

Existence and Uniqueness of Solutions to Ordinary Differential EquationsEdit

Peano’s theorem states that if   is continuous, then a solution to the ordinary differential equation   with initial condition   exists at least in some neighbourhood of time   – this solution need not be unique. Picard’s theorem states that if   is locally Lipschitz, then the solution for the ordinary differential equation   with initial condition   is unique when it exists. A comprehensive statement of these theorems is in Iserles[5], and there are proofs of these theorems in many books on ordinary differential equations (for example Birkhoff and Rota[6]).

Backward EulerEdit

We recall that the backward Euler method is given by   Note that if   is nonlinear, we need to solve a nonlinear equation in each step advancing the solution (numerical). It is usually hard to solve a nonlinear equation exactly using analytical methods, so we also use numerical methods. For our example equation, we get   This example has the advantage that we can find its solutions algebraically, so we can then examine the behavior of numerical schemes.

Convergence of Functional IterationEdit

We often use functional iteration to solve nonlinear equations. We recall that there are two popular methods: fixed-point iteration and Newton’s method.

Convergence of the Fixed-Point MethodEdit

We want to find a root of   We try to use the fixed-point method and to construct a sequence   where  .

Theorem Let   have a fixed-point  , be Lipschitz continuous for   with Lipschitz constant   and   be continuous on  . Then the fixed point method   converges to the unique fixed-point of   for  .

Proof Since   is Lipschitz continuous, we find that,   for  . Hence by induction we conclude that   Since  ,  , so we obtain a solution  , where   is the fixed point. We can show that the limit is unique by supposing that there are two different limits and reaching a contradiction. 

For a proof of the existence of the fixed-point under the assumptions used in this theorem, see a book on numerical analysis, such as Bradie[7] or Iserles[8].

Regarding our problem, we apply fixed-point iteration, we want to find the root of an equation of the form:

 

 

 

 

 

( 1)

When the timestep   is small enough then  . So fixed-point iteration is convergent provided the time-step is small enough. We note that eq. 1 has two roots, and so the domain of the initial iterate plays an important role in determining which root is choosen.

Convergence of Newton's MethodEdit

We now consider Newton's method. We want to find a root,   of   such that   Newton's method is a fixed-point method where the iterates are constructed by

 

 

 

 

 

( 2)

where  . If the function   is sufficiently well behaved, then Newton's method has a quadratic rate of convergence.

Theorem Suppose   is twice continuously differentiable and that its second derivative is bounded. Suppose also that there exists   for which  . Suppose   in the interval  ,   is finite in the same interval and   is small. Then, Newton's method is of quadratic convergence.

Proof   by Taylor expansion with Lagrange form remainder. In the above  . Since  , we have   so   Plug in the formula for  , from eq. 2 we have   Let   We have   and by our assumption, we know there is a constant   such that   Hence we have   for some finite constant  . So Newton's method is convergent provided   is sufficiently small. 

Regarding our problem, we consider  . Hence   and   is finite, so our problem satisfies all assumptions if we choose our initial data and initial iterates suitably. Hence the Newton iterations will converge and give an approximation to the nonlinear term in backward Euler's method.

Convergence of the Theta MethodEdit

The backward Euler, forward Euler and Crank-Nicolson methods are special case of the theta method, so we will first prove the convergence of the theta method to encompass these three methods. The theta method is the following algorithm,   where   and  . Notice that for   we obtain the Crank-Nicolson method or trapezoidal rule.

First, substituting the exact solution   and using the Taylor expansion we have

 

 

   

Subtracting the last expression from

 

we have that when   is small enough

  for  

  for  

where  . Using the triangle inequality and by the Lipschitz continuity of  , there exist constants   and   such that

  if  

  if  

When  , the theta method reduces to the trapezoidal rule. It is possible to show that the Crank-Nicolson method has second order convergence, see for example, Iserles[9]. Now let’s consider  ,  

We claim that  . We prove this statement by induction. When  ,  , since the initial conditions is exactly calculated. Now suppose this statement is true for  , where   and is a integer. We want to show this statement is true for  . Consider

 

then plug in

 

We have

 

 

So our claim is true for all   . Note that

 

by a Taylor expansion of the exponential function. Thus, we have

 

 

 

By our condition,   Therefore   So we have   and  . Hence the theta method is convergent of order 1 when  

Note that the backward Euler method is a special case of the theta method when  , so backward Euler's method is convergent of order 1. We arrive at our theorem.

Theorem Backward Euler’s method is convergent of order 1.

Remark If   is globally Lipschitz, then we can apply the above argument with respect to any time interval. If   is only locally Lipschitz, then we need to analyze the situation more carefully. First, by Picard’s theorem, there is a unique solution of this ordinary differential equation for a short amount of time. Indeed, we just need to know that the Lipschitz constant is finite without necessarily needing to know the exact value.

Remark If one did not know of Picard's theorem, one could deduce the existence and uniqueness of solutions to ODEs by using time discretization.

Now we consider   and  . The exact solution of this problem is  . So  . In our problem,   is clearly analytic and it is locally Lipschitz. It is easy to show   is not globally Lipschitz. If a function   is globally Lipschitz condition then there is a finite constant   such that   for all  . In our problem, let   and   it is easy to check  

We now discuss how one can find local Lipschitz constants  . When   is differentiable, we often just differentiate   and find the maximum value of its derivative in the domain of interest. In our example,   is simple and we only need to know that the Lipschitz constant is finite. So we use a more rough method to show that the Lipshitz constant is finite,

 

So it suffices to find the maximal value of   in this problem. In our problem,   is continuous. Furthermore,   will be positive all the time, because the initial value is positive and   is positive. A continuous function has finite maximal value in a closed and bounded set. Note that the exact solution of our problem is   , so  . So we know that the Lipschitz constant in our problem is finite.

Finally, we get the convergence of functional iteration and backward Euler’s method of our problem. Thus our numerical scheme for   with initial data   and   is convergent.

Corollary By the theorems for existence and uniqueness of the solution for ordinary differential equations and the Theorems on this page, we arrive at our final goal that the numerical solution generated by backward Euler's method with functional iteration exists and is unique when the time-step,   approaches zero.

Remark This requires careful choice of initial iterates when doing functional iteration.

Remark Typically, the exact solution of an ODE is not known, although it is possible to deduce local Lipschitz continuity. Should the solution become infinite, a numerical method will either not converge or display very large values if the approximate solution closely approximates the exact solution. Some care is required in interpreting such numerical simulations in these cases.

Example Programs which use Iteration to Solve a Nonlinear Ordinary Differential EquationEdit

The following two Matlab and Python programs demonstrate backward Euler’s method for the example equation. The first one uses fixed-point iteration to solve for the nonlinear term and the second one uses Newton’s method to solve for the nonlinear term.

 

 

 

 

( A)

A Matlab program to demonstrate fixed-point iteration Code download
% A program to solve y'=y^2 using the backward Euler
% method and fixed point iteration
% This is not optimized and is very simple

clear all; format compact; format short;
set(0,'defaultaxesfontsize',25,'defaultaxeslinewidth',.7,...
'defaultlinelinewidth',6,'defaultpatchlinewidth',3.7,...
'defaultaxesfontweight','bold')

n=10000; % Number of timesteps
Tmax=0.99; % Maximum time
y0=1; % Initial value
tol=0.1^10; % Tolerance for fixed point iterations
dt=Tmax/n; % Time step
y=zeros(1,n); % vector for discrete solution
t=zeros(1,n); % vector for times of discrete solution
y(1)=y0;
t(1)=0;
tic, % start timing
for i=1:n
    yold=y(i); ynew=y(i); err=1;
    while err>tol
        ynew=dt*yold^2+y(i);
        err=abs(ynew-yold);
        yold=ynew;
    end
    y(i+1)=ynew;
    t(i+1)=t(i)+dt;
end
toc, % stop timing
yexact=1./(1-t); max(abs(y-yexact)), % print the maximum error
figure(1); clf; plot(t,y,'r+',t,yexact,'b-.');
xlabel Time; ylabel Solution; legend('Backward Euler','Exact solution');
title('Numerical solution of dy/dt=y^2');

 

 

 

 

( Ap)

A Python program to demonstrate fixed-point iteration Code download
#!/usr/bin/env python
"""
A program to solve y'=y^2 using the backward Euler
method and fixed point iteration
This is not optimized and is very simple
"""

import time
import matplotlib.pyplot as plt

N = 1000	# Number of timesteps
tmax = 0.99	# Maximum time
y0 = 1
t0 = 0	# Initial value
tol = pow(0.1,10)	# Tolerance for fixed point iterations
h = tmax/N	# Time step

y = [y0]	# Variables holding the values of iterations
t = [t0]	# Times of discrete solutions



T0 = time.clock()
for i in xrange(N):
yold = y[i]
ynew = y[i]
err = 1
while err > tol:
ynew = h*pow(yold,2)+y[i]
err = abs(ynew-yold)
yold = ynew
y.append(ynew)
t.append(t[i]+h)

T = time.clock() - T0	
yexact = [1.0/(1.0-x) for x in t]

print
print "Exact value: y(%d)=%f" % (tmax, 1/(1-tmax))
print "Value given by aproximation: y(%d)=%f" % (tmax, y[-1])
maxerr=(max([abs(y[i] - yexact[i]) for i in xrange(len(y))]))
print "Maximum error: %f" % maxerr
print "Elapsed time is %f" % (T)

plt.figure()
plt.plot(t,y,'r+',t,yexact,'b-.')
plt.xlabel('Time')
plt.ylabel('Solution')
plt.legend(('Backward Euler', 'Exact solution'))
plt.title('Numerical solution of dy/dt=y^2')
plt.show()


 

 

 

 

( B)

A Matlab program to demonstrate Newton iteration Code download
% A program to solve y'=y^2 using the backward Euler
% method and Newton iteration
% This is not optimized and is very simple

set(0,'defaultaxesfontsize',25,'defaultaxeslinewidth',.7,...
'defaultlinelinewidth',6,'defaultpatchlinewidth',3.7,...
    'defaultaxesfontweight','bold')

n=100000; % Number of timesteps
Tmax=0.99; % Maximum time
y0=1; % Initial value
tol=0.1^10; % Tolerance for fixed point iterations
dt=Tmax/n; % Time step
y=zeros(1,n); % vector for discrete solution
t=zeros(1,n); % vector for times of discrete solution
y(1)=y0;
t(1)=0;
tic, % start timing
for i=1:n
    yold=y(i); ynew=y(i); err=1;
    while err>tol
        ynew=yold-(yold-y(i)-dt*yold^2)/(1-2*dt*yold);
        err=abs(ynew-yold);
        yold=ynew;
    end
    y(i+1)=ynew;
    t(i+1)=t(i)+dt;
end
toc, % stop timing
yexact=1./(1-t); max(abs(y-yexact)), % print maximum error
figure(1); clf; plot(t,y,'r+',t,yexact,'b-.');
xlabel Time; ylabel Solution;
legend('Backward Euler','Exact solution');
title('Numerical solution of dy/dt=y^2');

 

 

 

 

( Bp)

A Python program to demonstrate Newton iteration. Code download
#!/usr/bin/env python
"""
A program to solve y'=y^2 using the backward Euler
method and Newton iteration
This is not optimized and is very simple
"""

import time
import matplotlib.pyplot as plt

N = 1000 # Number of timesteps
tmax = 0.99 # Maximum value
t0 = 0 # Initial t value
y0 = 1 # Initial value y(t0) = y0
tol = 0.1**10 # Tolerance for fixed point iterations
h = (tmax - t0)/N # Time step

y = [y0] # List for discrete solutions
t = [t0] # List with grid of discrete values of t

T0 = time.clock() #Start timing

for i in xrange(N):
    yold = y[i]
    ynew = y[i]
    err =1
    while err > tol:
        ynew = yold-(yold-y[i]-h*(yold**2))/(1-2*h*yold)
        err = abs(ynew-yold)
        yold = ynew
    y.append(ynew)
    t.append(t[-1]+h)

T = time.clock() - T0 # Stop timing


print "Exact value y(%f) = %f " % (t[N], 1/(1-t[N]))
print "Value given by approx y_n(%f) = %f" % (t[N], y[N])
print "The error = y-y_n = %f " % (abs(1/(1-t[N]) - y[N]))
print "Time taken = %f " % (T)

yexact = [1.0/(1.0-x) for x in t]

plt.figure();
plt.plot(t,y,'r+',t,yexact,'b-.');
plt.xlabel('Time')
plt.ylabel('Solution')
plt.legend(('Backward Euler', 'Exact Solution'))
plt.title('Numerical solution of dy/dt=y^2')
plt.show()

ExercisesEdit

  1. Run the fixed-point iteration program in Matlab and check that the outcome is reasonable. Now investigate how changing the number of time steps taken to go from a time of 0 to a time of 0.99, and the tolerance for fixed point iterations affects the maximum error. In particular try a range of 1,000-1,000,000 (in powers of 10) for the number of time steps and a tolerance ranging from   (in powers of  ). You should observe that there is an ``ideal" combination of subdivisions and tolerance to minimize the error. What are these combinations? Do this whole process again using Newton iteration instead. How have the answers changed?
  2. Write a Matlab program to solve   with   using the Crank-Nicolson method and fixed point iteration. Explain why there are two fixed-points to which the fixed-point iteration can converge. Which of these fixed-points gives the correct approximation to the solution of the differential equation? Comment on how the choice of initial iterate for the fixed-point iteration determines the fixed-point to which the method converges.
    1. Show that the differential equation  , with   is not Lipschitz continuous.
    2. Find at least two analytical solutions to this differential equation.
    3. Compute a numerical solution to this differential equations using the forward Euler method.
    4. Compute a numerical solution to this differential equations using the backward Euler method. Be sure to try different initial guesses for the fixed-point iteration, not just the value at the previous time step; you should be able to calculate the influence of the choice of initial iterate on the selection of solution by the numerical method. Comment on this.
    5. Compute a numerical solution to this differential equations using the implicit midpoint rule. Be sure to try different initial guesses for the fixed point iteration, not just the value at the previous time step; you should be able to calculate the influence of the choice of initial iterate on the selection of “solution” by the numerical method. Comment on this.
    6. Repeat (d) and (e) with Newton iteration.
    7. Comment on the applicability of numerical methods for solving differential equations without unique solutions.
  1. Modify the program for the 1-D Allen-Cahn equation so that it uses the Crank-Nicolson and fixed-point iteration for the nonlinear term. You will need to calculate the nonlinear term in real space, so that your resulting scheme is   where   denotes the time step and   denotes the iterate. Stop the iterations once the maximum difference between successive iterates is sufficiently small.
  2. Modify the program for the 2-D Allen-Cahn equation so that it uses the Crank-Nicolson method and fixed-point iteration for the nonlinear term. You will need to calculate the nonlinear term in real space.

NotesEdit

  1. Iserles (2009)
  2. Lax, Burstein and Lax (1979)
  3. Hughes et al. (2008)
  4. Iserles (2009)
  5. Iserles (2009)
  6. Birkhoff and Rota (1989)
  7. Bradie (2006)
  8. Iserles (2009)
  9. Iserles (2009)

ReferencesEdit

Birkhoff, G; Rota, G.C. (1989). Ordinary Differential Equations (4 ed.). Wiley. 

Bradie, B. (2006). A Friendly Introduction to Numerical Analysis. Pearson. 

Hughes-Hallett, D.; Gleason A.M., Flath D.E., Lock P.F., Lomen D.O., Lovelock D., MacCallum W.G., Mumford D., Osgood B. G., Quinney D., Rhea K., Tecosky-Feldman J., Tucker T.W., and Bretscher O.K., Iovita A., Raskind W., Gordon S.P., Pasquale A., Thrash J.B. (2008). Calculus, Single and Multivariable (5th ed.). Wiley. 

Iserles, A. (2009). A First Course in the Numerical Analysis of Differential Equations. Cambridge University Press. 

Lax, A. (1976). Calculus with Applications and Computing. 1. Springer.