help using ODE45

cant figure out what im dong wrong

function [dxdt] = Alvarado_George_Dip(t,x,p,F)


m0=p(1) ;
m1= p(2) ;
m2 =p(3);
l1 =p(4);
l2=p(5);
g=p(6); 


x1=x(1);
x2=x(2);
x3=x(3);
x4=x(4);
x5=x(5);
x6=x(6);

temp = [((2*f-l1*(m1+2*m2)*(cos(x3)*dx4dt-sin(x3)*x4)-l2*m2*(cos(x5)*dx6dt-sin(x5)*x6))/(2*(m0+m1+m2)));...
    ((6*(g*(m1+2*m2)*sin(x3)-l2*m2*(cos(x3+x5)*dx6dt-sin(x3+x5)*x6^(2))+(m1+2*m2)*(x2*sin(x3)*x4-x2*sin(x3)*x4-cos(x3)*dx2dt)))/(3*l1*(m1+4*m2)+m2));...
    ((g*sin(x5)-l1*(cos(x3+x5)*dx4dt-sin(x3+x5)*x4^(2))-(x2+dx2dt)*sin(x5)*x6)/(2*(3*l2+1)))]  %[xddot;thetaddot1;thataddot2]
dx1dt = x2;
dx2dt = temp(1);
dx3dt = x4;
dx4dt = temp(2);
dx5dt=  x6
dx6dt= temp(3);

dxdt =[dx1dt;dx2dt;dx3dt;dx4dt;dx5dt;dx6dt];
end

then i use a call function

clear all;
close all;
clc;


m0=0.25 ;
m1= 0.1 ;
m2 =0.8;
l1 =0.25;
l2=0.2;
g=9.81; 
F=1;


p = [m0,m1,m2,l1,l2,g];
x0 = [0 0 0 0.01 0 0]';

%Calling Nonlinear Inverted Pendulum
[tout,xout] = ode45(@(t,x) Alvarado_George_Dip(t,x,p,F), [0 1], x0);

[James Tursa]

What is your specific question? Are you getting an error? (If so, please copy & paste the entire error message including the offending line). Are you getting unexpected results? (If so, please tell us why the results don't match your expectations).

[Stephan]

Can you provide the LaTex form of your system? If i correct the synax error in your code there still remains the problem, that your equations use variables before they are calculated.

Answer 1

[Alan Stevens]

Your temp equation can be written as the following three equations

dx2dt = k1*dx4dt + k2*dx6dt + k3
dx4dt = k4*dx2dt + k5*dx6dt + k6
dx6dt = k7*dx2dt + k8*dx4dt + k9

where the k's are the approprite combinations of x1, m1 , cos(x3), ... etc.

The equations can be expressed in matrix form as

M*dXdt = V; where
M = [1   -k1   -k2;
     -k4  1    -k5;
     -k7  -k8   1];
V = [k3; k6; k9];
dXdt = [dx2dt; dx4dt; dx6dt]

so you can solve these by

dXdt = M\V;
dx2dt = dXdt(1);
dx4dt = dXdt(2);
dx6dt = dXdt(3);

then you can proceed without calling any parameters before they are defined.

You will need to be very careful to ensure you define the k's correctly!

[George Alvarado]

[George Alvarado]

Maybe i approched the problem all wrong but what i want to is solve for x ,theta1, and theta 2

Answer 2

[Alan Stevens]

The equations above are a little different from the ones in your initial code, but I've used them in the coding below

m0=0.25 ;
m1= 0.1 ;
m2 =0.8;
l1 =0.25;
l2=0.2;
g=9.81; 
F=1;
p = [m0,m1,m2,l1,l2,g];
x0 = [0 0 0 0.01 0 0];

%Calling Nonlinear Inverted Pendulum
[tout,xout] = ode45(@(t,x) Alvarado_George_Dip(t,x,p,F), [0 0.5], x0);

theta0 = xout(:,1);
theta1 = xout(:,3);
theta2 = xout(:,5);
plot(tout,theta0,tout,theta1,tout,theta2),grid
xlabel('t'),ylabel('\theta')
legend('\theta_0','\theta_1','\theta_2')

function [dxdt] = Alvarado_George_Dip(~,x,p,F)

m0=p(1) ;
m1= p(2) ;
m2 =p(3);
l1 =p(4);
l2=p(5);
g=p(6); 

x1=x(1);
x2=x(2);
x3=x(3);
x4=x(4);
x5=x(5);
x6=x(6);

% Need to get dx2dt, dx4dt and dx6dt into separate equations 
% thet don't involve each other.
den2 = l1*(m1+2*m2)*cos(x3)+2*(m0+m1+m2);
den4 = 2*l1*(m1+3*m2);
den6 = l2*(m1+3*m2);
k2a = l2*m2*cos(x5)/den2;  % k2a*dx6dt
v2 = (2*F+l1*(m1+2*m2)*sin(x3)*x4^2+2*l2*m2*sin(x5)*x6^2)/den2;
k4a = 3*(m1+2*m2)*cos(x3)/den4;  % k4a*dx2dt
k4b = 3*l2*m2*cos(x3-x5)/den4;   % k4b*dx6dt
v4 = 3*(g*(m1+2*m2)*sin(x3)-l2*m2*sin(x3-x5)*x6^2)/den4;
k6a = 6*m2*cos(x5)/den6;       % k6a*dx2dt
k6b = 6*l1*m2*cos(x3-x5)/den6; % k6b*dx4dt
v6 = 6*m2*(g*sin(x5)+l1*sin(x3-x5)*x4^2)/den6;
%  1*dx2dt  +   0*dx4dt + k2a*dx6dt = v2
% k4a*dx2dt +   1*dx4dt + k4b*dx6dt = v4
% k6a*dx3dt + k6b*dx4dt +   1*dx6dt = v6
% The following three lines solve these
% to get dx2dt etc on their own
M = [1 0 k2a; k4a 1 k4b; k6a k6b 1];
V = [v2; v4; v6];
d246dt = linsolve(M,V);

dx1dt = x2;
dx2dt = d246dt(1);
dx3dt = x4;
dx4dt = d246dt(2);
dx5dt=  x6;
dx6dt= d246dt(3);

dxdt =[dx1dt;dx2dt;dx3dt;dx4dt;dx5dt;dx6dt];
end

Note that this blows up shortly after an end time of 0.5.

You will need to check my implementation carefully - I might have made a mistake in transferring the constants from the mathematical representation.