## Need help verifying code for Euler MEthod, RK2 and RK4

### David Geistlinger (view profile)

on 18 Nov 2019 at 23:29
Latest activity Edited by James Tursa

### James Tursa (view profile)

on 19 Nov 2019 at 1:56

### James Tursa (view profile)

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.

R2018b

### James Tursa (view profile)

on 19 Nov 2019 at 1:48
Edited by James Tursa

### James Tursa (view profile)

on 19 Nov 2019 at 1:56