How to solve s system of coupled nonlinear ODES already given the state vector

1 view (last 30 days)
I have a nonlinear dynamical system. I have forund the eqautions of motion, which results in the following state vector:
The variables I used are
I keep getting all kinds of errors. What am I foing wrong. Thank you!
function f = fun_1(t,x)
m0 = 0.25; m1 = 0.1; m2= 0.08;
L1 = 0.25; L2 = 0.2;
g = 9.81;
F = 0;
f=zeros(6,1);
f(1) =x(2);
f(2) =(-(L1*cos(x(3))*(m1+2*m2))/(2*(m0+m1+m2)))*f(4)+(-(m2*L2*cos(x(5)))/(2*(m0+m1+m2)))*f(6)+((L1*sin(x(3))*(m1+2*m2))/(2*(m0+m1+m2)))*x(4)^2+((m2*L2*sin(x(5)))/(2*(m0+m1+m2)))*x(6)^2+(2*F/(2*(m0+m1+m2)));
f(3) =x(4);
f(4) =(-(6*L1*cos(x(3))*(m1+2*m2))/(4*L1^2*(m1+3*m2)))*f(2)+(-(6*m2*L1*L2*cos(x(3)-x(5)))/(4*L1^2*(m1+3*m2)))*f(6)+(-(6*m2*L1*L2*sin(x(3)-x(5)))/(4*L1^2*(m1+3*m2)))*x(6)^2+((6*g*L1*(m1+2*m2)*sin(x(-(6*m2*L1*L2*sin(x(3)-x(5)))/3)))/(4*L1^2*(m1+3*m2)));
f(5) =x(6);
f(6) =(-(6*m2*L2*cos(x(5)))/(4*m2*L2^2))*f(2)+(-(6*m2*L1*L2*cos(x(3)-x(5)))/(4*m2*L2^2))*f(4)+(-(6*m2*L1*L2*sin(x(3)-x(5)))/(4*m2*L2^2))*x(4)^2+((6*m2*g*L2*sin(x(5)))/(4*m2*L2^2));
clc
clear all
close all
tspan = [0 10];
x0=[0;0;0;0;0;0];
[t, x] = ode23('fun_1', tspan, x0);
plot(t,x(:,3))
xlabel('t')
ylabel('theta(t)')

Accepted Answer

Alan Stevens
Alan Stevens on 23 Oct 2020
You need to solve for all xdots simultaneously - this is done in matrix format below. There were also a number of errors in your equations. The following works:
tspan = [0 10];
x0=[0;0;0;0;0;0];
[t, x] = ode23(@fun_1, tspan, x0);
plot(t,x(:,3))
xlabel('t')
ylabel('theta(t)')
function f = fun_1(t,x)
m0 = 0.25; m1 = 0.1; m2= 0.08;
L1 = 0.25; L2 = 0.2;
g = 9.81;
F = sin(t); % If F is zero and all the initial values of x are zero the result is
% inevitably a flat line! Replace with your own forcing function.
% Constants for use below
k24 = ((L1*cos(x(3))*(m1+2*m2))/(2*(m0+m1+m2)));
k26 = ((m2*L2*cos(x(5)))/(2*(m0+m1+m2)));
k42 = ((6*L1*cos(x(3))*(m1+2*m2))/(4*L1^2*(m1+3*m2)));
k46 = ((6*m2*L1*L2*cos(x(3)-x(5)))/(4*L1^2*(m1+3*m2)));
k62 = ((6*m2*L2*cos(x(5)))/(4*m2*L2^2));
k64 = ((6*m2*L1*L2*cos(x(3)-x(5)))/(4*m2*L2^2));
% LHS matrix of coefficients
M = [1 0 0 0 0 0;
0 1 0 k24 0 k26;
0 0 1 0 0 0;
0 k42 0 1 0 k46;
0 0 0 0 1 0;
0 k62 0 k64 0 1];
% RHS vector of constants
V = [x(2);
((L1*sin(x(3))*(m1+2*m2))/(2*(m0+m1+m2)))*x(4)^2+((m2*L2*sin(x(5)))/(2*(m0+m1+m2)))*x(6)^2+2*F/(2*(m0+m1+m2));
x(4);
-6*m2*L1*L2*sin(x(3)-x(5))/(4*L1^2*(m1+3*m2))*x(6)^2+6*g*L1*(m1+2*m2)*sin(x(3))/(4*L1^2*(m1+3*m2));
x(6);
6*m2*L1*L2*sin(x(3)-x(5))/(4*m2*L2^2)*x(4)^2 + 6*m2*g*L2*sin(x(5))/(4*m2*L2^2)];
% M*f = V where f is the vector of dxdt's
f = M\V;
end

More Answers (0)

Community Treasure Hunt

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

Start Hunting!