4.1 欧拉措施和龙格－库塔措施求微分方程数值解，画出解的图形，对成果进行比较分析如下方程供选择：d) 幂级数解解：原方程化为记，并用t替代x，则原方程化为：；且；于是可以用龙格－库塔法求解模型求解】：用Matlab作龙格－库塔法求解：% chapter 4—1.d%此函数是微分方程组function Xdot=ch41dfun(t,x)Xdot=[x(2),-cos(t)*x(1)]';%ch41d.mfunction I=ch41d(a)x0=[1,0]';[t,x]=ode45('ch41dfun',[0,a],x0);y=x(:,1);plot(t,y,'r');hold on;gtext('y');y1=1-1/jiecheng(2)*t.^2+2/jiecheng(4)*t.^4-9/jiecheng(6)*t.^6+55/jiecheng(8)*t.^8;plot(t,y1,'b');gtext('y1');[t,y,y1]hold off;运营程序可以得到 图1>> ch41d(10) 图2>> ch41d(15) 图3>> ch41d(20) 图4【成果分析】：由图1得：y(龙格－库塔措施)和y1(级数近似解)在0到大概1.5的区间内是完全吻合的，从x=1.5之后，两条曲线开始分离。
之后y的变化趋势可见图2－图4，呈幅度越来越大的上下震荡4.2 小型火箭初始重量为1 400kg，其中涉及1 080kg燃料火箭竖直向上发射时燃料的燃烧率为18kg/s，由此产生32 000N的推力，火箭引擎在燃料用尽时关闭设火箭上升时的空气阻力正比于速度的平方，比例系数为0.4kg/m求引擎关闭瞬间火箭的高度、速度、加速度，及火箭达到最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形解：解：先解出速度的微分方程如下：clearsyms a m t vm = 1400-18*t;a = (3-9.8*m-0.4*v^2)/m a =(18280+882/5*t-2/5*v^2)/(1400-18*t) 运用龙哥库塔算法画出从0时刻到60秒（燃料用尽）的图像odefun=@(t,v)(-((882*t)/5 - (2*v^2)/5 + 18280)/(18*t - 1400));[t,v]=ode45(odefun,[0:0.1:60],[0]);subplot(2,2,1)plot(t,v);grid ona=diff(v)/0.1;t2 = [0:0.1:59.9];subplot(2,2,2)plot(t2,a);grid ons = cumsum(v).*0.1;subplot(2,2,3)plot(t,s);grid on 因此v(60) = 267.3m/s a(60) = 0.9127m/s2 s(60) = 12200m4.3某容器盛满水后，底端直径为d0的小孔启动（如图1），根据水力学知识，当水面高度为h时，谁从小孔中流出的速度为v=0.6*（g*h）^0.5（其中g为重力加速度，0.6问哦小孔收缩系数）1） 若容器为倒圆锥形（如图1），现测得容器高和上底直径都为1.2m，小孔直径d为3cm，为水从小孔中流完需要多少时间；2min时水面高度是多少。
2） 若容器为倒葫芦形（如图2），现测得容器高1.2m，小孔直径d为3cm，由底端（记x=0）向上每隔0.1m测出容器的直径D（m）如表1所示，问水从小孔中流完需要多少时间；2min时水面高度是多少 图1X/m00.10.20.30.40.50.60.70.80.91.01.11.2D/m0.030.050.080.140.190.330.450.680.981.101.201.131.00 表1【模型分析】 由题知，水从小孔中流出，不仅与容器有关，还与水流速度v=0.6*（2*g*h）^0.5有关 第一小题容器是圆锥形，比较规则，但是由于水不断从小孔流出，容器中水的高度是不断变化的，水流速度没有一定的公式，因此要用到微积分解决由（1）知，水面直径等于水深水深为h时，流量为0.6(π/4)d^2*(gh)^0.5，0.6*(g*h)^(0.5)*π*(d0/2)^2*dt=π/4*h^2*dh则水深下降dh所需时间 ：dt=-[(π/4)h^2*dh]/[0.6(π/4)d^2*(gh)^0.5]=-[h^1.5*dh]/[0.6d^2*(g)^0.5]水深由1.2m至0定积分得水从小孔流完的时间：T（其中已知d=0.03m，g=9.8m*s(-2） 对于第二问：设两分钟（120S）后水深为X m ，由 dt=-[(π/4)h^2*dh]/[0.6*(π/4)*d^2*(gh)^0.5]=-[h^1.5*dh]/[0.6d^2*(g)^0.5] 则263.93-120 =X^2.5/[1.5*d^2*(g)^0.5] 以d=0.03m，g=9.8m*s(-2代入上式得 水深：X 第二小题容器为倒葫芦形，比较不规则，比较复杂，不仅要考虑水不断从小孔流出，容器中水的高度是不断变化的，水流速度没有一定的公式，因此要用到微积分解决，还要注意表1的倒葫芦形的不断变化，水深的高度变化是不规则的但仍可以用微积分。
由（2）知容器高1.2m，水深为h时，流量为0.6(π/4)d^2*(gh)^0.5，由于不同高度，倒葫芦形半径不同，用欧拉方程和龙格—库塔措施则水深下降dh所需时间 ：dt=t(k+1)-t(k)=-[(π/4)h^2*dh]/[0.6(π/4)d^2*(gh)^0.5]=-[h^1.5*dh]/[0.6d^2*(g)^0.5]然后运用循环for k=1：length(L)，t(k)=((h(k+1)-h(k))*(π/4)*d(k)^2)/(0.6*(π/4)*d^2*(g(1.2-h(k)))^0.5)，T=sum(t).可以求得水从小孔流完的总时间 对于第二问：设两分钟（120S）后水深为X m ，由 S=0,运用条件，当120-s<0.0001时s=s+t(k)，x(k) 把d=0.03m，g=9.8m*s(-2)代入上式得 水深：X 【模型求解】 开始条件h=1.2m，d0=1.2m，g=9.8m*s(-2)，d=0.03m 水深为h时，流量Q为0.6(π/4)d^2*(gh)^0.5，则水深下降dh所需时间 ：dt=-[(π/4)h^2*dh]/[0.6(π/4)d^2*(gh)^0.5]=-[h^1.5*dh]/[0.6d^2*(g)^0.5]水深由1.2m至0定积分得水从小孔流完的时间：T（其中已知d=0.03m，g-9.8m/s） 设两分钟（120S）后水深为X m ，由 dt=-[(π/4)h^2*dh]/[0.6*(π/4)*d^2*(gh)^0.5]=-[h^1.5*dh]/[0.6d^2*(g)^0.5] 则263.93-120 =X^2.5/[1.5*d^2*(g)^0.5] 以d=0.03，g=9.8代入上式得 水深：X 用欧拉方程和龙格—库塔措施则水深下降dh所需时间 ：dt=t(k+1)-t(k)=-[(π/4)h^2*dh]/[0.6(π/4)d^2*(gh)^0.5]=-[h^1.5*dh]/[0.6d^2*(g)^0.5] k1=0.15*sqrt(g*(x(n)))*d^2/(-43.6359*x(n)^8+213.0457*x(n)^7-414.873*x(n)^6+410.2075*x(n)^5-218.8936*x(n)^4+62.553*x(n)^3-8.3215*x(n)^2+0.49619*x(n)+.014892)^2; k2=0.15*sqrt(g*(x(n))-h*k1)*d^2/(-43.6359*(x(n)-h*k1)^8+213.0457*(x(n)-h*k1)^7-414.873*(x(n)-h*k1)^6+410.2075*(x(n)-h*k1)^5-218.8936*(x(n)-h*k1)^4+62.553*(x(n)-h*k1)^3-8.3215*(x(n)-h*k1)^2+0.49619*(x(n)-h*k1)+.014892)^2;x(n+1)=x(n)-h*(k1+k2)/2;【附录】(1)g=9.8;d=0.03;syms hy=-h^(1.5)/(0.6*d^2*sqrt(2*g));T=int(y,h,1.2,0);eval(T) x=((T-120)*(1.5*d^2*sqrt(2*g)))^(0.4);eval(x)运营成果如下>> sy1ans = 263.9316ans =0.9416(2)clear;g=9.8;d=0.03;k1=0;k2=0;h=0.4;x(1)=1.2;for n=1:1000 k1=0.15*sqrt(g*(x(n)))*d^2/(-43.6359*x(n)^8+213.0457*x(n)^7-414.873*x(n)^6+410.2075*x(n)^5-218.8936*x(n)^4+62.553*x(n)^3-8.3215*x(n)^2+0.49619*x(n)+.014892)^2; k2=0.15*sqrt(g*(x(n))-h*k1)*d^2/(-43.6359*(x(n)-h*k1)^8+213.0457*(x(n)-h*k1)^7-414.873*(x(n)-h*k1)^6+410.2075*(x(n)-h*k1)^5-218.8936*(x(n)-h*k1)^4+62.553*(x(n)-h*k1)^3-8.3215*(x(n)-h*k1)^2+0.49619*(x(n)-h*k1)+.014892)^2; x(n+1)=x(n)-h*(k1+k2)/2;end x(300) t=0:h:1000*h; plot(t,x); axis([0,400,0,1.21]);运营成果如下>> sy2ans =1.0278水从小孔流完需要263s，2min时水面高度是0.94m点头绪4.4 放射性废物的解决：将装满放射性废物的圆桶扔到水深300ft的海底，圆桶体积55gal，装满放射性废物时的圆桶重量为527.436lbf，在海中浮力为470.327。