how to convert ode15 code into Runge kutta 4th order

Question

Gur Pyari Rajput 2021-7-22

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/883719-how-to-convert-ode15-code-into-runge-kutta-4th-order

回答： prabhat kumar sharma 2024-6-4

CL= 1000;%pacing cycle length in ms

beats= 10;%number of beats in the simulation

options=[];%options for ode solver

for n=[1:beats]

[time X]=ode15s(@model,[0 CL],X0,options,1);

X0=X(size(X,1),:);

n %output beat number to the screen to monitor runtime progress

end

2 个评论
显示无隐藏无

James Tursa 2021-7-22

Ordinarily you could just replace ode15s with ode45, assuming your new problem isn't stiff. Or are you asking how to write your own RK4 code from scratch?

Gur Pyari Rajput 2021-7-29

Thanks James Tursa.... I will try...

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

prabhat kumar sharma 2024-6-4

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/883719-how-to-convert-ode15-code-into-runge-kutta-4th-order#answer_1467291

在 MATLAB Online 中打开

Hi Gur,

You can simply use ode45 instead of ode15 if your new problem is not stiff or If you are willing to write your own RK4 so you can follow the below refrence.

% Given parameters
CL = 1000; % Pacing cycle length in ms
beats = 10; % Number of beats in the simulation
X0 = [...]; % Initial conditions, define according to your model
% RK4 parameters
dt = 1; % Time step in ms, choose according to desired accuracy
% Placeholder for model function
% Ensure your model function is defined as: dxdt = model(t, X, ...)
% where dxdt is the rate of change of your variables (the derivative),
% t is the current time, and X is the current state vector
for n = 1:beats
    time = 0:dt:CL; % Time vector for the current beat
    X = zeros(length(time), length(X0)); % Initialize state matrix
    X(1, :) = X0; % Set initial state
    
    for i = 1:(length(time)-1)
        % Calculate the four slopes
        k1 = dt * model(time(i), X(i, :)', 1);
        k2 = dt * model(time(i) + dt/2, X(i, :)' + k1/2, 1);
        k3 = dt * model(time(i) + dt/2, X(i, :)' + k2/2, 1);
        k4 = dt * model(time(i) + dt, X(i, :)' + k3, 1);
        
        % Update the state
        X(i+1, :) = X(i, :) + (k1 + 2*k2 + 2*k3 + k4)'/6;
    end
    
    X0 = X(end, :); % Prepare initial condition for the next beat
    disp(['Beat number: ', num2str(n)]); % Output beat number to monitor progress
end