Need help Euler Method

I need help with my code for programming this second order Euler equation. I tried to insert new fuction into the euler code I already had and not sure when im going wrong.

Here is the question:

Here is my code:

% The system of figure shown below consists of a uniform rigid link of
% total mass m and length L, two linear springs of stiffnesses k1 and k2,
% and a viscous damper of damping constant c.  When the link is horizontal,
% the springs are unstretched.  The motion of the system is described by
% the following second order differential equation:
%
% 1/3mL^2theta2 + k1L^2(1-cos(theta))sin(theta) + k2L^2sin(theta)cos(theta)
% - mgL/2cos(theta) + cL^2(theta1)cos^2(theta) = 0 

% 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

h = 0.005;
N = 5;
theta(1) = pi/10;
theta1(1) = 0;

for n=1:N
    x(n+1) = n*h;
    theta(n+1) = theta(n) + h*(1/3*m*L^2*theta2 + k1*L^2*(1-cos(theta))*sin(theta) + k2*L^2*sin(theta)*cos(theta) - m*g*L/2*cos(theta) + c*L^2*(theta1)*(cos(theta))^2) + theta(n);
end
plot(x,theta)
format long
disp (theta)

Not sure what I am doing wrong. Any help would be great. Thank you in advance

Answer 1

[James Tursa]

The fundamental thing you are doing wrong is that you don't have the proper size state vector. For a 1st order equation, the state vector has one element. For a 2nd order equation, the state vector has two elements. If you had two coupled 2nd order equations, the state vector has 2x2 = 4 elements. Etc.

In your case, you have a single 2nd order equation, so your state vector will have two elements. I will keep your same variable names for this discussion. You are using theta(1), theta(2), etc. as your solution states:

theta(1) is the solution at time x(1)
theta(2) is the solution at time x(2)
etc.

But this is inadequate as you only have a single element state at each time point. You need two elements at each time point (one for theta and one for theta_dot) because of your 2nd order equation. I find it convenient to use a matrix to hold the solution with the columns as the states at the time points (but you could easily use rows if you wanted). In this case the states would be:

theta = zeros(2,N); % Preallocate the result
theta(1:2,1) is the solution at time x(1) for theta and theta_dot
theta(1:2,2) is the solution at time x(2) for theta and theta_dot
etc.

I.e., theta(1,1) is theta at time x(1) and theta(2,1) is theta_dot at time x(1)

So your initialization equations of:

theta(1) = pi/10;
theta1(1) = 0;

would instead be:

theta(1,1) = pi/10;
theta(2,1) = 0;

And your state update equation of:

theta(n+1) = theta(n) + h*(1/3*m*L^2*theta2 + k1*L^2*(1-cos(theta))*sin(theta) + k2*L^2*sin(theta)*cos(theta) - m*g*L/2*cos(theta) + c*L^2*(theta1)*(cos(theta))^2) + theta(n);

would instead be something like:

theta(:,n+1) = theta(:,n) + h*(stuff);

The "stuff" would be a 2x1 vector. The first element of this vector will be theta_dot at the current time, and the second element will be theta_doubledot at the current time. This is close to what you have already, but don't use plain theta without subscripts in the right hand side expression since theta is the ENTIRE solution matrix! Instead, use an expression based on the current solution state which is theta(:,n). (for you to figure out)

And to plot just the "theta" angle part and not the "theta_dot" part, do this:

plot(x,theta(1,:))

Also, the following is not correct:

N = 5;

The 5 value is the end time, not the number of time points. So this should be something like this instead:

tstop = 5; % end time
N = round(tstop/h) + 1;

Give all this a shot and come back when you have more problems with your updated code.