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:piz0=[ 0; 0; 10*sin(i); 10*cos(i)];fun=@(t,z)[z(3);z(4);..

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:piz0=[ 0; 0; 10*sin(i); 10*cos(i)];fun=@(t,z)[z(3);z(4);..
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内,有多少条“线”

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:piz0=[ 0; 0; 10*sin(i); 10*cos(i)];fun=@(t,z)[z(3);z(4);..
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可以减少计算量