Need help verifying code for Euler MEthod, RK2 and RK4

9 views (last 30 days)
This is the problem I am trying to solve:
Project 1.PNG
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.

Accepted Answer

James Tursa
James Tursa on 19 Nov 2019
Edited: James Tursa on 19 Nov 2019
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.

More Answers (0)

Categories

Find more on Programming in Help Center and File Exchange

Products


Release

R2018b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!