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

12 views (last 30 days)
David Geistlinger on 15 Nov 2019
Answered: Thomas Satterly on 15 Nov 2019
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.

Thomas Satterly on 15 Nov 2019
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)