Contents
Download: MATLAB files
See deriv.m
Euler method
clear all xi = 0.5; A = [0 1; -1 -2.0*xi]; B = [0; 1]; N = 21; u = linspace(0,10,N); x = [0;0]; du = u(2)-u(1); y = zeros(size(u)); for n=1:N y(n) = x(1); dx = (A*x + B)*du; x = x + dx; end plot(u,y);
Exact solution (underdamped xi <1)
xi = 0.5; vf = 1; wn = sqrt(1-xi^2); A = [ 1 0; -xi wn ]; b = [-vf; 0]; K = A\b; ye = vf + (K(1)*cos(wn*u) + K(2)*sin(wn*u)).*exp(-xi*u); plot(u,ye);
Runge-Kutta
x = [0;0]; yr = zeros(size(u)); err = zeros(size(u)); h = u(2) - u(1); for n=1:N t = u(n); yr(n) = x(1); k1 = deriv(t, x) * h; k2 = h * deriv(t + h/3.0, x + k1/3.0); k3 = h * deriv(t + h/3.0, x + k1/6.0 + k2/6.0); k4 = h * deriv(t + h/2.0, x + k1/8.0 + (3.0/8.0)*k3); k5 = h * deriv(t + h, x + k1/2.0 - 1.5 * k3 + 2.0 * k4); errx = (2.0*k1 - 9.0*k3 + 8.0*k4 - k5)/30.0; err(n) = errx(1); dx = (k1 + 4.0*k4 + k5)/6.0; x = x + dx; end plot(u,yr);
MATLAB ODE45
ysolve = ode45(@deriv,[0 10],[0;0]); tm = ysolve.x; ym = ysolve.y(1,:); plot(tm,ym);
Compare exact with Euler
subplot(2,1,1); plot(u, y, u, ye); xlabel('\omega t'); ylabel('response'); subplot(2,1,2); plot(u,y-ye); xlabel('\omega t'); ylabel('error');
Compare exact with Runge-Kutta
subplot(2,1,1); plot(u, yr, u, ye); xlabel('\omega t'); ylabel('response'); subplot(2,1,2); plot(u,yr-ye,u,err); xlabel('\omega t'); ylabel('error');