How to solve s system of coupled nonlinear ODES already given the state vector
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Missael Hernandez
le 22 Oct 2020
Réponse apportée : Alan Stevens
le 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 commentaires
Réponse acceptée
Alan Stevens
le 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 commentaires
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!