MATH 128A, SUMMER 2009: PROGRAMMING ASSIGNMENT 3: SOLUTIONS
(1) Let’s change the notation for Newton’s method so it doesn’t conflict with the notation for ODEs:
Newton’s method finds a root of g(x) using the iteration pj+1 =pj−g(pj)
g0(pj).
The equation for wi+1 is equivalent to g(x) = 0, where g(x) = wi+h
2[f(ti+1, x) + f(ti, wi)]−x.
To use Newton’s method, we need this function’s derivative: g0(x) = h
2
∂f
∂y (ti+1, x)−1. Write down:
p0=wi
pj+1 =pj−wi+h
2[f(ti+1, pj) + f(ti, wi)] −pj
h
2
∂f
∂y (ti+1, pj)−1
(2) In implicit_trap.m:
function [ti,wi]=implicit_trap(f,dfdy,a,b,y0,N,tol,max_iterations)
h = (b-a)/N;
ti = linspace(a, b, N+1);
wi(1) = y0;
for i=1:N
% Set up a root-finding problem, g(w(i+1)) = 0.
g = @(x) -x + wi(i) + h/2 * ( f (ti(i+1),x) + f(ti(i),wi(i)));
dg = @(x) -1 + h/2 * (dfdy(ti(i+1),x) );
wi(i+1) = newton(g, dg, wi(i), tol, max_iterations);
end
(3) f = @(t,y) y.*(10-y); dfdy = @(t,y) 10-2*y; y0 = 2;
tExact = linspace(0, 4);
wExact = 10 ./ (1 + 4*exp(-10*tExact));
[tEuler, wEuler] = euler(f, 0, 4, y0, 4*4);
[tITrap, wITrap] = implicit_trap(f, dfdy, 0, 4, y0, 4*4, 1e-4, 20);
plot(tExact, wExact, tEuler, wEuler, tITrap, wITrap)
legend(’Exact’, ’Euler’, ’ITrap’,0)
Figure 1. Plot for problem 3 (MATLAB)
0 0.5 1 1.5 2 2.5 3 3.5 4
2
3
4
5
6
7
8
9
10
11
12
Exact
Euler
ITrap
Date: Wednesday 8/12.
1
