Trapezoid algorithm on an ODE
6 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Drake
am 17 Feb. 2021
Kommentiert: Drake
am 17 Feb. 2021
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.
0 Kommentare
Akzeptierte Antwort
James Tursa
am 17 Feb. 2021
Bearbeitet: James Tursa
am 17 Feb. 2021
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
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!