Solving a System of 2nd Order Nonlinear ODEs

Hi, I need help with setting up the code for the below 2 non-linear differential equations. I do not have symbolic toolbox. Below is how I tried it but didn't work. Any help will be appreciated.

clear all;clc
function dydt=mbd(M,m,g,lc,k)

y1=y(1);
y2=y(2);
z1=z(1); %z1
z2=z(2); %z2
y5=(1/(M+m))*((m*l/2)*(sin(z(1)))*y6 + (m*l/2)*(cos(z(1)))*(z(2))^2 ...
               -c*y(2)-k*y(1));
y6=(3/(m*l^2))*((m*l/2)*(sin(z(1)))*y5 - (m*g*l/2)*(sin(z(1))));

dydt=[y1;y2;y3;y4;y5;y6];
end

Answer 1

[Sam Chak]

Hey @mpz

It doesn't work that way because you have created an algebraic loop error where depends y5 that is your , which depends on y6 that is from the beginning.

I'll try to explain how to solve it without using difficult math or matrix operations. Here, the governing equations are given by

Equation 5:

Equation 6:

.

Equation 6 can be rearranged to become

so that it can substituted into Equation 5 to eliminate

and then be simplified to become

... Equation (7).

Note the the right-hand side of Equation 7 does not have the term . Subtituting Equation 7 into the simplified Equation 6 yields

... Equation (8).

Now you can rewrite the mbd ode function in state-space form correctly. Note that the ODEs should have only 4 states (instead of 6 states):

If you find my technical explanation is helpful, please vote and accept the answer.

[mpz]

@Sam Chak Thanks. I wrote the function and call to solve it but getting errors. Is my syntax incorrect? As you can see I have 4 initial conditions (x,xdot,theta,thetadot). I'm used to 2 initial conditions so I'm not used to 4.

% function
function dydt=mbd(y,M,m,g,l,C,k)
dydt=zeros(6,1);
dydt(1)=y(1); %linear displacement
dydt(2)=y(2); %linear speed
dydt(3)=y(3); %angular displacement
dydt(4)=y(4); %angular speed
dydt(5)=(1/(M+m-(3/4)*(m*(sin(y(3)))^2)))*(((m*l/2)*(cos(y(3))*(y(4))^2)) ...
     - ((3/4)*(m*g*((sin(y(3)))^2))) - (C*(y(2))) - (k*(y(1))));
dydt(6)=((3*sin(y(3)))/(2*l*(M+m-(3/4)*(m*(sin(y(3)))^2))))*(((m*l/2)*(cos(y(3))*(y(4))^2)) ...
    - ((3/4)*(m*g*((sin(y(3)))^2))) - (C*(y(2))) - (k*(y(1))))- ((3*g*sin(y(3)))/(2*l));
% dydt=[y1;y2;y3;y4;y5;y6];
end

% Call m file
clear all;clc;close all
g=0.2;
M=2;
m=pi()/4;
l=1;
k=0.5;
C=0.25;
t=0:0.001:20;
initial_x = 0.5 ;
initial_dxdt = 0 ;
initial_theta = pi()/3;
initial_dthetadt = 0;
[t,y] = ode45(@mbd,t,[initial_x initial_dxdt initial_theta initial_dthetadt]);
p=y(:,1); %% displacement (rad)
q=y(:,2); %% angular speed (rad/s)
figure (1)
plot (t,p)

Error

Index exceeds the number of array elements (1).

Error in mbd (line 5)
dydt(2)=y(2); %linear speed

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in mbdc (line 13)
[t,y] = ode45(@mbd,t,[initial_x initial_dxdt initial_theta initial_dthetadt]);
 
>> 
 

[Torsten]

You have 4 equations and unknowns. Why do you address 6 in "mbd" ?

[mpz]

@Torsten I see what you mean but I don't know how do to setup the xdotdot[dydt(5)] and thetadotdot[dydt(6)]. Do I setup it as follows?

function dydt=mbd(y,M,m,g,l,C,k)
dydt=zeros(6,1);
dydt(1)=y(2); %linear speed
dydt(2)=y(4); %angular speed
dydt(3)=(1/(M+m-(3/4)*(m*(sin(y(3)))^2)))*(((m*l/2)*(cos(y(3))*(y(4))^2)) ...
     - ((3/4)*(m*g*((sin(y(3)))^2))) - (C*(y(2))) - (k*(y(1))));
dydt(4)=((3*sin(y(3)))/(2*l*(M+m-(3/4)*(m*(sin(y(3)))^2))))*(((m*l/2)*(cos(y(3))*(y(4))^2)) ...
    - ((3/4)*(m*g*((sin(y(3)))^2))) - (C*(y(2))) - (k*(y(1))))- ((3*g*sin(y(3)))/(2*l));
% y(1) linear displacement and y(3) angular displacement
end

[Torsten]

function dydt=mbd(y,M,m,g,l,C,k)
% y1 = x
% y2 = theta
% y3 = xdot
% y4 = thetadot
dydt=zeros(4,1);
dydt(1)=y(3); %linear speed
dydt(2)=y(4); %angular speed
dydt(3)=(1/(M+m-(3/4)*(m*(sin(y(2)))^2)))*(((m*l/2)*(cos(y(2))*(y(4))^2)) ...
    - ((3/4)*(m*g*((sin(y(2)))^2))) - (C*(y(3))) - (k*(y(1))));
dydt(4)=((3*sin(y(2)))/(2*l*(M+m-(3/4)*(m*(sin(y(2)))^2))))*(((m*l/2)*(cos(y(2))*(y(4))^2)) ...
    - ((3/4)*(m*g*((sin(y(2)))^2))) - (C*(y(3))) - (k*(y(1))))- ((3*g*sin(y(2)))/(2*l));
end

[mpz]

@Torsten thanks. When I used my call code with that function, I got the error below.

clear all;clc;close all
g=0.2;
M=2;
m=pi()/4;
l=1;
k=0.5;
C=0.25;
t=0:0.001:20;
initial_x = 0.5 ;
initial_dxdt = 0 ;
initial_theta = pi()/3;
initial_dthetadt = 0;
[t,y] = ode45(@mbd,t,[initial_x initial_dxdt initial_theta initial_dthetadt]);
p=y(:,1); %% displacement (rad)
q=y(:,2); %% angular speed (rad/s)
figure (1)
plot (t,p)

Error

Index exceeds the number of array elements (1).

Error in mbd (line 8)
dydt(1)=y(3); %linear speed

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in mbdc (line 13)
[t,y] = ode45(@mbd,t,[initial_x initial_dxdt initial_theta initial_dthetadt]);
 
>> 

[Torsten]

[t,y] = ode45(@(t,y)mbd(y,M,m,g,l,C,k),t,[initial_x initial_dxdt initial_theta initial_dthetadt]);

instead of

[t,y] = ode45(@mbd,t,[initial_x initial_dxdt initial_theta initial_dthetadt]);

How should MATLAB know what additional parameters you want to transfer to mbd if you don't specify them ?

[mpz]

Thank you

[Sam Chak]

Hi @mpz

Sorry, a little late. The code should look like this. Also thank @Torsten for showing the corrections.

function mpz

close all

clc

tspan = [0 20]; % time span of simulation

y0 = [0.5 pi/3 0 0]; % initial values

[t, y] = ode45(@(t, y) odefcn(t, y), tspan, y0);

plot(t, y, 'linewidth', 1.5)

grid on

xlabel('Time, t [sec]')

ylabel({'$x,\; \theta,\; \dot{x},\; \dot{\theta}$'}, 'Interpreter', 'latex')

title('Time responses of the system states')

legend({'x', '$\theta$', '$\dot{x}$', '$\dot{\theta}$'}, 'Interpreter', 'latex', 'location', 'best')

end

function dydt = odefcn(t, y) % Ordinary Differenctial Equations

dydt = zeros(4,1);

g = 0.2;

M = 2;

m = pi/4;

l = 1;

C = 0.25;

k = 0.5;

dydt(1) = y(3);

dydt(2) = y(4);

dydt(3) = (1/(M + m - (3/4)*(m*(sin(y(2)))^2)))*(((m*l/2)*(cos(y(2))*(y(4))^2)) - ((3/4)*(m*g*((sin(y(2)))^2))) - (C*(y(3))) - (k*(y(1))));

dydt(4) = ((3*sin(y(2)))/(2*l*(M + m - (3/4)*(m*(sin(y(2)))^2))))*(((m*l/2)*(cos(y(2))*(y(4))^2)) - ((3/4)*(m*g*((sin(y(2)))^2))) - (C*(y(3))) - (k*(y(1)))) - ((3*g*sin(y(2)))/(2*l));

end

Result:

[mpz]

Thank you @Sam Chak

Answer 2

[Alan Stevens]

You haven't passed y and z to the funcion in

function dydt=mbd(M,m,g,lc,k)

Probably needs to be more like

function dydt=mbd(t,y,z,M,m,g,lc,k)

with the calling function modified accordingly.