Numerical Solutions of Differential Equations
Demos 1:Euler
% euler numerical integration of y' = cos(x).
% here's an inline function... saves an mfile
% for easy little expressions...
dsin = 'cos(x)';
dy = inline(dsin, 'x'); % all args are strings
y = zeros(1,100); % answers to be computed
N = 100; % number of steps
InitVal = 0; % between 0
FinalVal = 2*pi; % and 2pi
Range = FinalVal-InitVal;
h = Range/N; %step size
x = 0.0;
y(1) = 0; %initial condition
for i = 2:100
x = x+h;
y(i) = y(i-1)+ h*dy(x);
end
plot(y, 'r');
% or sort of vectorizing, get something that's
% hard to turn into experimental code
xy = linspace(0,2*pi); % full of the x(i) first
for i = 2:100
xy(i) = xy(i-1)+ cos(xy(i)); % replaced by y(i)
end
figure
plot(y, 'b')
Next, Euler code for the mass-spring-damper system showing the
integration of a second-order system expressed as two coupled
first-order
equations.
% euler numerical integration of mx'' -kx' +sx = g(t)
function msby = eulmsb(m,k,s)
N = 500 % steps and answers
y1 = zeros(1,N); % y(t) answers to be computed
y2 = zeros(1,N); % intermediate y'(t) answers computed
InitVal = 0; % between 0
FinalVal = 100; % and here
Range = FinalVal-InitVal;
h = Range/N; %step size
x = 0.0;
y1(1) = 10; % initial condition for y
y2(1) = 0; % initial condition for y'
y2a = 0; % initial condition for y' not saving y'(t)
for i = 2:N
x = x+h;
% y1(i) = y1(i-1)+ h*y2(i-1);
% y2(i) = y2(i-1)+ h*(y1(i-1)*(-s/m) +y2(i-1)*(k/m)) +g(i-1);
y2a = y2a+ h*(y1(i-1)*(-s/m) +y2a *(k/m)) +g(i-1);
y1(i) = y1(i-1)+ h*y2a;
end
plot(y1, 'r');
msby = y1;
end
function driven = g(k)
driven= sin((k/10)*(2*pi));
% driven = 0;
end
Function Handles
abs(-4) % call abs function, returns 4
@abs % a fn. handle. returns @abs
my_abs = @abs % now I have name for handle
my_abs(-5) % returns 5 (? maybe shouldn't work?)
@abs(-5) % syntax (? semantics surely? ) error
[T,Y] = ode23(@abs,timespan,y0) % OK
[T,Y] = ode23(my_abs,timespan,y0) % OK, same thing
% this next implements the lambda expression above
function apply (funh,x,y)
funh(x,y)
end
%so
apply(@max,5,4) % returns 5
% but....(!)
@+ % syntax error (??-- OK in Lisp)
Matlab ODE23
For the Lorenz attractor.
function dy = lorenz(t, y)
%Lorenz Attractor. t a scalar, y a vector of y's
%dy returned as a column vector of dy's
% y(1)... are the x,y,z
dy = zeros(3,1);
dy(1) = 10*(y(2) - y(1));
dy(2) = - y(1)*y(3) + 28*y(1) - y(2);
dy(3) = y(1)*y(2)-(8/3)*y(3);
end
% Script to Run Lorenz
lorenzhandle = @lorenz;
y0 = [10, -10, 20]; % initial conditions
timespan = linspace(0, 20, 1500); %specify 1500 instants
[T,Y] = ode23 (lorenzhandle,timespan, y0); % boom
plot(Y(:,1),Y(:,3), 'r');
hold on
y0 = [10, -10, 20.1]; % initial conditions close to previous
[T,Y] = ode23 (lorenzhandle,timespan, y0);
plot(Y(:,1),Y(:,3), 'b'); %BUT big difference in the two solns
Back to Slides
Last update: 04/20/2011: RN