清空筛选条件
清空筛选条件

Trapezoid algorithm on an ODE

6 次查看（过去 30 天）
Drake
Drake 2021-2-17
评论： Drake 2021-2-17
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
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 个）

类别

在 Help Center 和 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