我正在嘗試使用 MATLAB 中的一些函式來數值求解一對耦合的二階常微分方程
\ddot{x} = f(x,y,\dot{x},\dot{y})
\ddot{y} = f(x,y,\dot{x},\dot{y})。
我能夠讓它僅與一個二階 ODE 一起作業，但我嘗試使用的代碼不適用于一對 ODE。
函式 odeToVectorField 有效地采用二階 ODE 并將其寫為一對耦合一階 ODE 的向量。ode45 是通常的 Runge-Kutta 求解方法。xInit 和 yInit 對應于 x 和 y 的初始條件，然后目標是在特定時間間隔內繪制 x 和 y 與時間的關系圖。
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)
uj5u.com熱心網友回復：
只是為了比較，不使用符號運算式，可以將這個等式實作為
function dV = M(t,V)
gamma1=0.1;
gamma2=0.1;
a=1;
m=1;
g=9.8;
d=1;
x = V(1); dx = V(2); y = V(3); dy = V(4);
ddx = (gamma1*dx)/(a m*d^2 (m/2)*d^2*cos(y-x)) (gamma2*dy)/(a (m/2)*cos(y-x)) - ( (m/2)*d^2*sin(y-x)*(dx^2 - dy^2))/(a m*d^2 (m/2)*d^2*cos(y-x)) - ((m/2)*d^2*dx^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));
ddy = (gamma1*dx)/((m/2)*d^2*cos(y-x)) (gamma2*dy)/a - ( (m/2)*d^2*sin(y-x)*(dx^2 - dy^2))/((m/2)*d^2*cos(y-x)) - ((m/2)*d^2*dx^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;
dV = [dx ddx dy ddy];
end%function
interval = [0 20];
xInit = [2 0];
yInit = [2 0];
vSol = ode45(M,interval,[ xInit yInit]);
tValues = linspace(0,20,100);
xValues = deval(vSol,tValues,1);
plot(tValues,xValues)
這有效，但報告t=0.244和之間的奇點t=0.245。
轉載請註明出處，本文鏈接：https://www.uj5u.com/yidong/389987.html