Need help with heuns method (rk2) for second order DE

Question

David Geistlinger 2019-11-15

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/491200-need-help-with-heuns-method-rk2-for-second-order-de

回答： Thomas Satterly 2019-11-15

采纳的回答： Thomas Satterly

在 MATLAB Online 中打开

Here is the probelm:

Here is my Matlab Code:

% Runge-Kutta Second Order Method
clc
clear
% Parameters
g = 9.81;
m = 3;
L = 1;
k1 = 100;
k2 = 150;
c = 1.5;
% Inputs
h = 0.005;
t_final = 5;
N = t_final/h;
% Initial Conditions
t(1) = 0;
x(0) = pi/10;
v(0) = 0;
% Update loop
for i=1:N
    t(i+1) = t(i)+h;
    xs=x(i)+h*(v(i));
    vs=v(i)+h*(((-3/m)*((k1*(1-cos(x(i))*sin(x(i))) + (k2*sin(x(i))*cos(x(i))) + (c*v*cos(x(i))^2)))) + ((3*g*cos(x(i)))/(2*L)));
    x(i+1)=x(i)+h/2*(v(i)+vs);
    v(i+1)=v(i)+h/2*((((-3/m)*((k1*(1-cos(x(i))*sin(x(i))) + (k2*sin(x(i))*cos(x(i))) + (c*v*cos(x(i))^2)))) + ((3*g*cos(x(i)))/(2*L)))-((-3/m)*((k1*(1-cos(xs)*sin(xs)) + (k2*sin(xs)*cos(xs)) + (c*v*cos(xs)^2)))) + ((3*g*cos(xs))/(2*L)));
end
figure(1); clf(1)
plot(t,x)

Getting a few errors:

1) Unable to perform assignment because the left and right sides have a different number of elements.

Error in ProjectRk2 (line 43)

x(i+1)=x(i)+h/2*(v(i)+vs);

2) Array indices must be positive integers or logical values.

Error in ProjectRk2 (line 33)

x(0) = pi/10;

Any help with this code would be nice or if there is a easier way to do rk2 and rk4 i would be interested in hearing. My knowledge isnt the best with matlab but not the worst.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Thomas Satterly 2019-11-15

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/491200-need-help-with-heuns-method-rk2-for-second-order-de#answer_401733

在 MATLAB Online 中打开

1) You forgot to index the variable "v" in a couple places

2) Matlab indicies start at 1, not 0.

Corrected code:

% Runge-Kutta Second Order Method
clc
clear
% Parameters
g = 9.81;
m = 3;
L = 1;
k1 = 100;
k2 = 150;
c = 1.5;
% Inputs
h = 0.005;
t_final = 5;
N = t_final/h;
% Initial Conditions
t(1) = 0;
x(1) = pi/10;
v(1) = 0;
% Update loop
for i=1:N
    t(i+1) = t(i)+h;
    xs=x(i)+h*(v(i));
    vs=v(i)+h*(((-3/m)*((k1*(1-cos(x(i))*sin(x(i))) + (k2*sin(x(i))*cos(x(i))) + (c*v(i)*cos(x(i))^2)))) + ((3*g*cos(x(i)))/(2*L)));
    x(i+1)=x(i)+h/2*(v(i)+vs);
    v(i+1)=v(i)+h/2*((((-3/m)*((k1*(1-cos(x(i))*sin(x(i))) + (k2*sin(x(i))*cos(x(i))) + (c*v(i)*cos(x(i))^2)))) + ((3*g*cos(x(i)))/(2*L)))-((-3/m)*((k1*(1-cos(xs)*sin(xs)) + (k2*sin(xs)*cos(xs)) + (c*v(i)*cos(xs)^2)))) + ((3*g*cos(xs))/(2*L)));
end
figure(1); clf(1)
plot(t,x)