MATLABで4つの結合微分方程式を解く

Question

MATLABで解決しようとしている結合ODEのセットがあります。方程式を以下に示す。

I 4つの境界条件を有する：X（0）、Y（0）、V（0）、シータ（0）。 dsolveでこれを解決しようとすると、明確な解決策が見つからないという警告が表示されます。 ここに私が使ったコードがあります。

syms x(t) y(t) v(t) theta(t) 

% params 
m = 80;    %kg 
g = 9.81;   %m/s^2 
c = 0.72;   % 
s = 0.5;   %m^2 
theta0 = pi/8;  %rad 
y0 = 0;    %m 
rho = 0.94;   %kg/m^3 


% component velocities 
xd = diff(x,t) == v*cos(theta); 
yd = diff(y,t) == v*sin(theta); 

% Drag component 
D = c*rho*s/2*(xd^2+yd^2); 

% Acceleration 
vd = diff(v,t) == -D/m-g*sin(theta); 

% Angular velocity 
thetad = diff(theta,t) == -g/v*cos(theta); 

cond = [v(0) == 10,y(0) == 0, x(0) == 0, theta(0) == theta0]; 
dsolve([xd yd vd thetad],cond) 

Answer 1

これは少しいくつかの並べ替えの振り子のように見える...

あなたの方程式は、少なくとも、振り子のODEに似ている用語

dθ(t)/dt = C·cos(θ(t)) 

を持って、それが同じ問題を抱えています。この方程式の閉じた形の解は知られていません。私はそれが存在しないことが証明されていると信じていますが、私は100％確実ではありません。

とにかく、数値的にはケーキです。ここode45でそれを行う方法の例を示します（変数tspanでキャプチャ）厳密解とは異なり、あなたは開始時刻と終了時刻を指定する必要がありますことを

function my_ode() 

    % parameters 
    m = 80;  % kg 
    g = 9.81; % m/s² 
    c = 0.72; % - 
    s = 0.5; % m² 
    rho = 0.94; % kg/m³ 

    theta0 = pi/8; % rad 
    v0  = 10; % m/s 
    x0  = 0; % m 
    y0  = 0; % m 

    tspan = [0 10]; % s 


    % function to compute derivative of 
    % Z = [x, y, th, v] 
    function dZdt = odefcn(~,Z) 

     % rename for clarity (NOTE: x and y are not used) 
     th = Z(3); cth = cos(th); 
     v = Z(4); sth = sin(th); 

     % Compute derivatives 
     dxdt = v * cth; 
     dydt = v * sth; 
     dthdt = -g/v * cth;   
     dvdt = -c*rho*s*v^2/(2*m) - g*sth; 

     % assign to ouptut respecting either row or columnvector inputs 
     dZdt = Z; 
     dZdt(:) = [dxdt dydt dthdt dvdt];   

    end 

    % Integrate the ODE 
    Z0 = [x0 y0 theta0 v0];  
    [t,Z] = ode45(@odefcn, tspan, Z0); 



    % Example outputs 
    x = Z(:,1); th = Z(:,3); 
    y = Z(:,2); v = Z(:,4); 

    F = figure; hold on 
    P = get(F, 'position');  
    set(F, 'position', [P(1:2) 3*P(3) P(4)]);  

    subplot(1,3,1) 
    plot(x,y, 'b'), grid on 
    xlabel('x [m]'), ylabel('y [m]') 

    subplot(1,3,2) 
    plot(t,v, 'b'), grid on 
    xlabel('t'), ylabel('v [m/s]') 

    subplot(1,3,3) 
    plot(t,th, 'b'), grid on 
    xlabel('t'), ylabel('\theta [rad]') 

end 

注意。 Dを簡略化するために、ID cos²θ + sin²θ = 1を使用したことにも注意してください。

出力例：