Решите 4 связанных дифференциальных уравнения в MATLAB

Question

У меня есть набор связанных ODE, которые я хочу решить с помощью MATLAB. Уравнения приведены ниже.

У меня 4 граничные условия: х (0), у (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

Это выглядит немного как маятник какой-то ...

Ваше уравнение имеет термин

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

, который похож на ОДУ маятника, по крайней мере, его имеет ту же проблему: решение замкнутой формы для этого уравнения неизвестно. Я считаю, что даже доказано, что этого не существует, но я не уверен на 100%.

Во всяком случае, численно это кусок торта. Вот пример того, как сделать это с ode45:

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 

Обратите внимание, что в отличие от точного решения, вы должны будете указать время начала и окончания (захваченное в переменных tspan). Обратите также внимание на то, что я использовал идентификатор cos²θ + sin²θ = 1 для упрощения D.

Пример вывода: