Need help verifying code for Euler MEthod, RK2 and RK4

This is the problem I am trying to solve:

Here is my Euler Method Code:

% Euler Method

clc
clear 
format long

% Assume:
g = 9.81; % m/s^2
m = 3; % kg
L = 1; % m
k1 = 100; % N/m
k2 = 150; % N/m
c = 1.5; % N*s/m

f=@(t,y) [y(2);-(3*(2*L*k1*sin(y(1)) - g*m*cos(y(1)) + 2*L*c*cos(y(1))^2*y(2) - 2*L*k1*cos(y(1))*sin(y(1)) + 2*L*k2*cos(y(1))*sin(y(1))))/(2*L*m)];
disp('Solution is')
[T,Y]=eulerSystem(f,[0,5],[pi/10;0],0.005);
plot(T,Y(1,:));
title('Plot of th(t)');
figure;
plot(T,Y(2,:));
title('Plot of th''(t)');
function [t,y]=eulerSystem(Func,Tspan,Y0,h)
t0=Tspan(1);
tf=Tspan(2);
N=(tf-t0)/h;
y=zeros(length(Y0),N+1);
y(:,1)=Y0;
t=t0:h:tf;
for i=1:N
y(:,i+1)=y(:,i)+h*Func(t(i),y(:,i));
end
end

Here is code for Rk2:

% 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)

Here is code for Rk4:

% Runge-Kutta Fourth Order Method
clc
clear
% Parameters
g = 9.81;
m = 3;
L = 1;
k1 = 100;
k2 = 150;
c = 1.5;
dt = 0.005;
t = 0:dt:5;
x0 = pi/10;
v0 = 0;

% Rk4
f1 = @(t,x,v) v;
f2 = @(t,x,v) (((-3/m)*((k1*(1-cos(x)*sin(x)) + (k2*sin(x)*cos(x)) + (c*v*cos(x)^2)))) + ((3*g*cos(x))/(2*L)));
h = dt;

x = zeros(length(t),1);
v = zeros(length(t),1);
x(1) = x0;
v(1) = v0;

for i = 1:length(t)-1
    K1x = f1(t(i),x(i),v(i));
    K1v = f2(t(i),x(i),v(i));
    K2x = f1(t(i) + h/2, x(i) + h*K1x/2, v(i) + h*K1v/2);
    K2v = f2(t(i) + h/2, x(i) + h*K1x/2, v(i) + h*K1v/2);
    K3x = f1(t(i) + h/2, x(i) + h*K2x/2, v(i) + h*K2v/2);
    K3v = f2(t(i) + h/2, x(i) + h*K2x/2, v(i) + h*K2v/2);
    K4x = f1(t(i) + h, x(i) + h*K3x, v(i) + h*K3v);
    K4v = f2(t(i) + h, x(i) + h*K3x, v(i) + h*K3v);
    x(i+1) = x(i) + h/6*(K1x + 2*K2x + 2*K3x + K4x);
    v(i+1) = v(i) + h/6*(K1v + 2*K2v + 2*K3v + K4v);
end

plot(t,x)

All codes work, however getting very different figures which I think should be similiar. I have to use these methods and not ODE. I am pretty sure the Euler method code is right and the RK2 code looks similiar but the Rk4 code is does not show a figure I can justify. Please help any advice with the codes especcially the Rk4 one is most appreicated.

Answer 1

[James Tursa]

Your biggest problem is that you don't have the differential equations coded correctly. And the main cause of that is because you picked a bad way to code them. The first thing you did was great, which was to make a function handle f for the Euler case that returned the entire state derivative, a 2x1 vector. You coded that correctly, and things plotted well.

But then ...

You abandoned that correctly coded function handle and started writing the derivative code again from scratch multiple times. And, it turns out, you coded it incorrectly this time.

You can go over those derivative lines one by one and double check the coding on each one (turns out the (1-cos(theta))*sin(theta) term you have coded incorrectly as (1-cos(theta)*sin(theta))), or you can take my suggestion and scrap all of that code and instead call the same function handle in your RK2 and RK4 cases that you used in your Euler case. This should be the way it was coded from the beginning, having your derivative code all in one place so that each method is guaranteed to be using the same derivative. And if you change your RK2 and RK4 code to use vectors for your states instead of individual x and v states, the coding will be easier and you will be set up to easily use ode45 as a double check on your results.

Finally, it was difficult to compare your coded derivative with the assignment text because you changed things around and the equation was lengthy. IMO, the equations were complicated enough to warrant a separate multi-line function rather than a function handle. It would have made it much easier to debug.