![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/834035/image.png)
Numerically solving a pair of coupled second order ODES with odeToVectorField
6 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Hollis Williams
le 14 Déc 2021
Commenté : Nikhil Kapoor
le 16 Avr 2022
I am attempting to use some of the functions in MATLAB to numerically solve a pair of coupled second order ODEs of the form
\ddot{x} = f(x,y,\dot{x},\dot{y})
\ddot{y} = f(x,y,\dot{x},\dot{y}).
I am able to get it to work with just one second-order ODE, but the code I am trying to does not work for a pair of ODEs.
The function odeToVectorField effectively takes a second order ODE and writes it as a vector for a pair of coupled first order ODEs. ode45 is the usual Runge-Kutta solution method. xInit and yInit correspond to the initial conditions for x and y and the aim is then to plot both x and y against time and also x against y over time.
gamma1=0.1;
gamma2=0.1;
a=1;
m=1;
g=9.8;
d=1;
syms x(t) y(t)
eqn1=diff(x,2)== (gamma1*diff(x))/(a + m*d^2 + (m/2)*d^2*cos(y-x)) + (gamma2*diff(y))/(a+ (m/2)*cos(y-x)) - ( (m/2)*d^2*sin(y-x)*(diff(x)^2 - diff(y)^2))/(a + m*d^2 + (m/2)*d^2*cos(y-x)) - ((m/2)*d^2*diff(x)^2*(y-x))/(a+ (m/2)*cos(y-x)) - ((m/2)*d*(3*g*sin(x) + g*sin(y)))/(a + m*d^2 + (m/2)*d^2*cos(y-x)) - ((m/2)*d*g*sin(y))/(a+ (m/2)*cos(y-x))
eqn2=diff(y,2)== (gamma1*diff(x))/((m/2)*d^2*cos(y-x)) + (gamma2*diff(y))/a - ( (m/2)*d^2*sin(y-x)*(diff(x)^2 - diff(y)^2))/((m/2)*d^2*cos(y-x)) - ((m/2)*d^2*diff(x)^2*(y-x))/a - ((m/2)*d*(3*g*sin(x) + g*sin(y)))/((m/2)*d^2*cos(y-x)) - ((m/2)*d*g*sin(y))/a
V = odeToVectorField(eqn1,eqn2)
M = matlabFunction(V,'vars',{'t','Y'})
interval = [0 20];
xInit = [2 0];
yInit = [2 0];
ySol = ode45(M,interval,xInit, yInit);
tValues = linspace(0,20,100);
yValues = deval(ySol,tValues,1);
plot(tValues,yValues)
0 commentaires
Réponse acceptée
Star Strider
le 14 Déc 2021
‘... the code I am trying to does not work for a pair of ODEs.’
Yes, it does. Look at the ‘Subs’ result from odeToVectorField to see that everything is there.
The problem that remains is that this is now a
degree system so ‘xInit, yInit’ need to be concatenated with square brackets for it to work —
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/834035/image.png)
gamma1=0.1;
gamma2=0.1;
a=1;
m=1;
g=9.8;
d=1;
sympref('AbbreviateOutput',false);
syms x(t) y(t)
eqn1=diff(x,2)== (gamma1*diff(x))/(a + m*d^2 + (m/2)*d^2*cos(y-x)) + (gamma2*diff(y))/(a+ (m/2)*cos(y-x)) - ( (m/2)*d^2*sin(y-x)*(diff(x)^2 - diff(y)^2))/(a + m*d^2 + (m/2)*d^2*cos(y-x)) - ((m/2)*d^2*diff(x)^2*(y-x))/(a+ (m/2)*cos(y-x)) - ((m/2)*d*(3*g*sin(x) + g*sin(y)))/(a + m*d^2 + (m/2)*d^2*cos(y-x)) - ((m/2)*d*g*sin(y))/(a+ (m/2)*cos(y-x))
eqn2=diff(y,2)== (gamma1*diff(x))/((m/2)*d^2*cos(y-x)) + (gamma2*diff(y))/a - ( (m/2)*d^2*sin(y-x)*(diff(x)^2 - diff(y)^2))/((m/2)*d^2*cos(y-x)) - ((m/2)*d^2*diff(x)^2*(y-x))/a - ((m/2)*d*(3*g*sin(x) + g*sin(y)))/((m/2)*d^2*cos(y-x)) - ((m/2)*d*g*sin(y))/a
[V,Subs] = odeToVectorField(eqn1,eqn2)
M = matlabFunction(V,'vars',{'t','Y'})
interval = [0 20];
xInit = [2 0];
yInit = [2 0];
ySol = ode45(M,interval,[xInit, yInit]);
figure
plot(ySol.x, ySol.y)
grid
legend(string(Subs), 'Location','best')
tValues = linspace(0,20,100);
yValues = deval(ySol,tValues,1);
plot(tValues,yValues)
There are still problems, however the code essentially works.
.
1 commentaire
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!