clc<br>clear all<br>f=dsolve('cos(t)*(sin(pi*t/T)^2)==D2y');<br>syms T theta R y0 y0d x0 x0d t x y<br>s=diff(f,t);<br>f0=subs(f,T,35);<br>syms C5 C4<br>C5=solve(subs(f0,t,0)==y0,C5);<br>df=diff(f,t);<br>dff=subs(df,T,35);<br>C4=solve(subs(dff,t,0)==y0d,C4);<br>f1=eval(f);<br>f2=x0+t*x0d;<br>f3=subs(f1,t,(x-x0)/x0d);<br>s0=subs(eval(s),t,(x-x0)/x0d);<br>u=subs(s0,T,35);<br>u=subs(u,y0,R*sin(theta));<br>u=subs(u,x0,R*cos(theta));<br>u=subs(u,R,1);<br>f4=subs(f3,y0,R*sin(theta));<br>f4=subs(f4,x0,R*cos(theta));<br>f4=subs(f4,T,35);<br>f4=subs(f4,R,1);<br>d=1000;<br>c=0;<br>h=1;<br>for i=1:100<br>x0d=2*(rand-0.5);<br>y0d=2*(rand-0.5);<br>theta=rand*pi;<br>g=eval(f4);<br>u0=eval(u);<br>if atan(abs(cot(theta)))+atan(abs(y0d/x0d))<=pi/2<br>c(h)=solve(g==d);<br>m(h)=subs(eval(u),x,c(h));<br>XD(h)=x0d;<br>h=h+1;<br>end<br>end<br>cmin=min(c);<br>cmax=max(c);<br>k=linspace(-10000,10000,1000);<br>cc=hist(c,k);<br>cc=cc/length(c);<br>bar(k,cc)<br>dd=hist3(XD,m,k);<br>bar3(XD,m) 您好,由于内容包含色情低俗信息,建议亲尽量不要发布这种信息哦~ 把我暂存的matlab码hx了。。。