Asked by David Geistlinger
on 18 Nov 2019 at 23:29

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 by James Tursa
on 19 Nov 2019 at 1:48

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

Accepted Answer

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.

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.