How to solve s system of coupled nonlinear ODES already given the state vector
1 view (last 30 days)
Show older comments
Missael Hernandez
on 22 Oct 2020
Answered: Alan Stevens
on 23 Oct 2020
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)')
0 Comments
Accepted Answer
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
0 Comments
More Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!