Trapezoid algorithm on an ODE

t(1,1) = 0;    % Starting time (arbitrary units)
u(1,1) = 1;    % Initial position (arbitrary units)
u(1,2) = 0;    % Initial velocity (arbitrary units)
tend = 4*pi; % Simulation run time.
Nsteps = 200;  % Number of integration steps;
dt = (tend-t(1))/Nsteps;  % Time step size (arbitrary units)

%% Set up the differential equation
% This is the "right hand side" function of the differential equation, T'
k = 1; m = 1; omega = sqrt(k/m);       % Physical parameters of oscillator
dudt = @(u) [ u(2) -omega^2 * u(1) ];  % Right-hand side function

%% Observe for a specified amount of time
for idx = 1:tend/dt                     
    t(idx+1,1) = t(idx,1) + dt;        
    utmp = u(idx,:) + dt*dudt(u(idx,1));
    u(idx+1,:) = u(idx,:) + dt*(dudt(t(idx,:))+dudt(utmp))/2;
end

I am getting an error with this code that does an ODE with a trapezoid method, any help will be appreciated.

Answer 1

[James Tursa]

You need to pass the entire state to the derivative function, not just one element of the state. E.g.,

utmp = u(idx,:) + dt*dudt(u(idx,1));

should be

utmp = u(idx,:) + dt*dudt(u(idx,:));

and you need to replace your t(idx,:) with u(idx,:). So this

u(idx+1,:) = u(idx,:) + dt*(dudt(t(idx,:))+dudt(utmp))/2;

should be

u(idx+1,:) = u(idx,:) + dt*(dudt(u(idx,:))+dudt(utmp))/2;

To make things clear that you want the function handle to return a 1x2 vector, maybe include the comma separator:

dudt = @(u) [ u(2) , -omega^2 * u(1) ];  % Right-hand side function

[Drake]

Thank you so much!