Question

MATLAB问题求解：微分方程组解的概率分布

m=4.1888*10^(-6);k=9.3*10^(-7);g=9.8;%参数
ts=0:0.001:2;%时间区间
%初始值 x(0)=0 y(0)=0 x'(0)=10*sin(i) y(0)'=cos(i)
for i=0:0.025:pi
z0=[ 0; 0; 10*sin(i); 10*cos(i)];
fun=@(t,z)[z(3);z(4);...
-k/m*z(3)*sqrt(z(3)^2+z(4)^2);...
-k/m*z(4)*sqrt(z(3)^2+z(4)^2)-g];
[t,z]=ode45(fun,ts,z0);
subplot(121),plot(t,z(:,1:2));
xlabel('t');legend('x','y');
axis([0 1.48 0 7 ])
hold on
subplot(122),plot(z(:,1),z(:,2));
axis([-5 5 0 3 ])
xlabel('x');ylabel('y');
hold on
end

我解出这个方程组之后，图像大家可以看看。
想知道当y=0时，x的密度函数。就是x在每个小区间x±δt内，有多少条“线”

Answer 1

function []=test2

m=4.1888*10^(-6);k=9.3*10^(-7);g=9.8;%参数

N=5000;  %角度分段

ts=[0 2];%时间区间，确保时间足够回到y=0就可以

fun=@(t,z)[z(3);z(4);...

    -k/m*z(3)*sqrt(z(3)^2+z(4)^2);...

    -k/m*z(4)*sqrt(z(3)^2+z(4)^2)-g];

op=odeset('Events',@evenfun);

ang=linspace(0,pi,N);

x=nan(size(ang));

for i=1:length(ang);

    %初始值 x(0)=0 y(0)=0 x'(0)=10*sin(i) y(0)'=cos(i)

    z0=[    0;      0;   10*sin(ang(i)); 10*cos(ang(i))];

    [t,z,TE,ZE]=ode45(fun,ts,z0,op);%TE触发事件时间，ZE事件被触发时的状态

    %当检查到y重新回到0，停止计算，这时候TE是发生的时刻，ZE是这时候的状态

    %可以ZE(2)在误差范围内为0表示回到y=0，ZE(1)就是这时候的横坐标

    %t=0时刻y=0，为了排除一开始的点，下边设定条件要时间TE大于0的值才能算回到y=0

    %平抛和向下抛的过程都不可能回到y=0了

    if TE>eps,x(i)=ZE(1);end

end

x=x(~isnan(x));%得到的向量x就是每条轨迹除了起始点外和y=0的交点的x值

%分区间统计，可以统计不同区间内交点的个数

hist(x,50);%这里给出分50个区间统计的结果

%角度分段N越大统计越精确，不过计算需要时间也长

%统计区间越多得到的分布分辨越好，但是质量越差

%通常根据N的取值，确定分区间数，太多太少都不好

%最后图像是统计频数的分布

end

function [value,isterminal,direction]=evenfun(t,z)

    value=z(2);%事件z(2)=0,也就是y=0

    isterminal=1;%事件触发，停止计算

    direction=-1;%方向负，也就是从y大于0到y小于0

end

这里将角度分了5000分，可能需要计算半分钟左右的时间

如果将所有的统计加起来应该只有一半也就是2500的事件

因为从平抛开始后边的都不可能再回到y=0

所以实际上，可以将角度范围控制在0~pi/2可以减少计算量