小程序：matlab解微分方程

function f = ball_trajectory(t, y);


g = 9.81;
pi = 3.14159265;
air_density = 1.20;


M = 57.5 * 0.001;
D = 64.0 * 0.001;
R = D * 0.5;
A = pi * R * R;


Cd = 0.507;
ratio_Cl_to_S = 0.58;


F = 0.5 * Cd * air_density * A / M;
B = 0.5 * ratio_Cl_to_S * air_density * A * R / M;


V = sqrt(y(4) * y(4) + y(5) * y(5) + y(6) * y(6));
W = sqrt(y(7) * y(7) + y(8) * y(8) + y(9) * y(9));


f = zeros(9, 1);


f(1) = y(4);
f(2) = y(5);
f(3) = y(6);


f(4) = -1.0 * F * V * y(4) + B * (y(8) * y(6) - y(9) * y(5));
f(5) = -1.0 * F * V * y(5) + B * (y(9) * y(4) - y(7) * y(6));
f(6) = -1.0 * F * V * y(6) + B * (y(7) * y(5) - y(8) * y(4)) - g;


%f(7) = 0;
%f(8) = 0;
%f(9) = 0;


%f(7) = -10.0 * F * V * y(7);
%f(8) = -10.0 * F * V * y(8);
%f(9) = -10.0 * F * V * y(9);


f(7) = -0.2 * y(7);       % 0.2
f(8) = -0.2 * y(8);       % 0.2
f(9) = -0.2 * y(9);       % 0.2


%f(7) = -0.2;      
%f(8) = -0.2;

%f(9) = -0.2;

%f(7) = 0;
%f(8) = 0;
%f(9) = 0;


%f(7) = -10.0 * F * V * y(7);
%f(8) = -10.0 * F * V * y(8);
%f(9) = -10.0 * F * V * y(9);


% dw/dt = 2.9 * 10^(-4) / R * V * w


f(7) = -2.9 * 10^(-4) * V * y(7);
f(8) = -2.9 * 10^(-4) * V * y(8);
f(9) = -2.9 * 10^(-4) * V * y(9);


%f(7) = -0.2;      
%f(8) = -0.2;
%f(9) = -0.2;