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')

%  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

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