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.
采纳的回答
James Tursa
2021-2-17
编辑：James Tursa
2021-2-17
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
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
