Trapezoid algorithm on an ODE
33 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponse acceptée
James Tursa
le 17 Fév 2021
Modifié(e) : James Tursa
le 17 Fév 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
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!