UCB Math 128A, Spring 2009: Programming Assignment 4
Due May 6
Consider the initial-value problem
y0=f(t, y), a ≤t≤b, y(a) = α. (1)
The Backward Euler method is given by
w0=α(2)
wi+1 =wi+hf(ti+1, wi+1) (3)
Since fdepends on the unknown wi+1, (3) is an equation that must be solved at each time step, and
in general fmight be a non-linear function. In this assignment, we will implement the backward
Euler method in MATLAB using Newton’s method for solving (3). The solver will then be used to
solve a stiff differential equation.
1. Write down Newton’s method for solving (3) for wi+1, using the initial guess w(0)
i+1 =wi. The
iterations will depend on the partial derivative ∂f/∂y =fy(t, y ).
2. Implement a MATLAB function backeuler.m of the form
function [t,w] = backeuler(f,dfdy,a,b,alpha,N,maxiter,tol)
which implements the backward Euler method using the Newton’s method from 1. Here, fis
the function f(t, y), dfdy is the function fy(t, y), Nis the number of time steps, maxiter is the
maximum number of iterations in the Newton solver (the code should break if this is exceeded
without convergence), and tol is the tolerance between two successive Newton iterations. At
each time step, print the value of the Newton updates to monitor the convergence (similar to
the newton table.m function on the web page).
Use the following simple test to check your solver:
f=@(t,y) -y*sin(y);
df=@(t,y) -sin(y)-y*cos(y);
a=0; b=3; alpha=1; N=20; maxiter=20; tol=1e-12;
[t1,w1]=backeuler(f,df,a,b,alpha,N,maxiter,tol);
[t2,w2]=rk4(f,a,b,alpha,N);
plot(t1,w1,t2,w2)
legend(’Backward Euler’,’Runge-Kutta 4’)
The two curves should be close to each other, and if the number of steps Nis increased they
should get even closer (the RK4 method is much more accurate and can be considered the
true solution). Do not include this test in your report, it is only meant as a verification of
your backward Newton method.
Turn page −→
