Download Math 151B Homework 4: Comparing Numerical Methods for Solving Differential Equations and more Assignments Mathematics in PDF only on Docsity!
Math 151B Homework #4 – not turned in
- Adams-Bashforth two point method.
Solve the initial value problem
y′^ = −ty, 0 ≤ t ≤ 4 , y(0) = 1
using the Adams-Bashforth two point method.
Do the same with the midpoint method.
Answer:
function [t,w] = ab2(fname,a,b,alpha,h) f = str2func(fname); N = (b-a)/h; t(1) = a; w(1) = alpha;
% find a second point using midpoint method w(2) = w(1) + 0.5hf(t(1),w(1)); w(2) = w(1) + h*f(t(1)+(h/2),w(2)); t(2) = t(1) + h;
% now do the rest with AB- for i = 2:N w(i+1) = w(i) + 0.5h( 3*f(t(i),w(i)) - f(t(i-1),w(i-1)) ); t(i+1) = t(i) + h; endfor endfunction
0
1
0 0.5 1 1.5 2 2.5 3 3.5 4
Exact Euler Midpt AB 2pt
Both of these are 2nd^ order methods. Which one requires more function evaluations?
Answer: As written, both methods use two function evaluations per step. However, the AB-2 method could store up the value of f (t(i), w(i)) for use on the next step, thus removing a function call.
- Second order equation
The motion of a pendulum of length ` is described by the equation
θ′′^ = −
g `
sin θ
where θ is the angular deviation from vertical.
Let g/` = 1, so we have θ′′^ = − sin θ.
a. Linearize this equation (i.e. approximate sin θ with a linear function) and solve exactly.
Answer: θ′′^ = −θ which has general solution θ(t) = a 1 sin(t) + a 2 cos(t). With initial angle θ 0 and zero initial velocity, we get θ(t) = θ 0 sin(t).
d. Challenge! Write an implicit Euler’s method to solve this system.
Answer:
function [t,w1,w2] = implicit_euler_2d(a,b,alpha,beta,h) N = (b-a)/h; tol = 0.0001; t(1) = a; w1(1) = alpha; w2(1) = beta; for i = 1:N t(i+1) = t(i) + h; w1i = w1(i); w2i = w2(i); w1t = w1i; w2t = w2i; for j=1: numer = w1t - w1i - h(w2i-hsin(w1t)); denom = 1 - (h^2)cos(w1t); w1(i+1) = w1t - numer/denom; if abs(w1(i+1)-w1t)<tol break; else w1t = w1(i+1); endif if j> printf("can’t meet tolerance\n"); return; endif endfor for j=1: numer = w2t - w2i + hsin(w1i+hw2t); denom = 1 - (h^2)cos(w1i+h*w2t); w2(i+1) = w2t - numer/denom; if abs(w2(i+1)-w2t)<tol break; else w2t = w2(i+1); endif if j> printf("can’t meet tolerance\n"); return; endif endfor endfor endfunction
