數(shù)學實驗課程設計常微分方程數(shù)值解

1、 數(shù)學實驗報告1. 題目： 某容器盛滿水后，底端直徑為d0的小孔開啟（如圖1），根據水力學知識，當水面高度為h時，誰從小孔中流出的速度為v=0.6*（g*h）0.5（其中g為重力加速度，0.6問哦小孔收縮系數(shù)）1） 若容器為倒圓錐形（如圖1），現(xiàn)測得容器高和上底直徑都為1.2m，小孔直徑d為3cm，為水從小孔中流完需要多少時間；2min時水面高度是多少。2） 若容器為倒葫蘆形（如圖2），現(xiàn)測得容器高1.2m，小孔直徑d為3cm，由底端（記x=0）向上每隔0.1m測出容器的直徑D（m）如表1所示，問水從小孔中流完需要多少時間；2min時水面高度是多少。 圖1X/m00.10.20.30.40.5

2、0.60.70.80.91.01.11.2D/m0.030.050.080.140.190.330.450.680.981.101.201.131.00 表12. 分析： 由題知，水從小孔中流出，不僅與容器有關，還與水流速度v=0.6*（2*g*h）0.5有關。 第一小題容器是圓錐形，比較規(guī)則，但是由于水不斷從小孔流出，容器中水的高度是不斷變化的，水流速度沒有一定的公式，所以要用到微積分解決。由（1）知，水面直徑等于水深。水深為h時，流量為0.6(/4)d2*(gh)0.5，0.6*(g*h)(0.5)*(d0/2)2*dt=/4*h2*dh則水深下降dh所需時間 ：dt=-(/4)h2*dh

3、/0.6(/4)d2*(gh)0.5=-h1.5*dh/0.6d2*(g)0.5水深由1.2m至0定積分得水從小孔流完的時間：T（其中已知d=0.03m，g=9.8m*s(-2） 對于第二問：設兩分鐘（120S）后水深為X m ，由 dt=-(/4)h2*dh/0.6*(/4)*d2*(gh)0.5=-h1.5*dh/0.6d2*(g)0.5 則263.93-120 =X2.5/1.5*d2*(g)0.5 以d=0.03m，g=9.8m*s(-2代入上式得 水深：X 第二小題容器為倒葫蘆形，比較不規(guī)則，比較復雜，不僅要考慮水不斷從小孔流出，容器中水的高度是不斷變化的，水流速度沒有一定的公式，所

4、以要用到微積分解決，還要注意表1的倒葫蘆形的不斷變化，水深的高度變化是不規(guī)則的但仍可以用微積分。 由（2）知容器高1.2m，水深為h時，流量為0.6(/4)d2*(gh)0.5，由于不同高度，倒葫蘆形半徑不同，用歐拉方程和龍格庫塔方法則水深下降dh所需時間 ：dt=t(k+1)-t(k)=-(/4)h2*dh/0.6(/4)d2*(gh)0.5=-h1.5*dh/0.6d2*(g)0.5然后利用循環(huán)for k=1：length(L)，t(k)=(h(k+1)-h(k)*(/4)*d(k)2)/(0.6*(/4)*d2*(g(1.2-h(k)0.5)，T=sum(t).可以求得水從小孔流完的總時

5、間。 對于第二問：設兩分鐘（120S）后水深為X m ，由 S=0,利用條件，當120-s<0.0001時s=s+t(k)，x(k) 把d=0.03m，g=9.8m*s(-2)代入上式得 水深：X 相關資料：龍格庫塔方法求解龍格庫塔方法思想：用v,v上若干個點的倒數(shù)，對他們做線性組合得到平均斜率，就可能得到更高階的精度，這就是龍格庫塔方法思想。為了演算方便本課程設計采用二階龍格庫塔公式，即求v,v上的二個倒數(shù)，將他們加權平均得平均斜率。我們設有形如： 其中函數(shù)滿足HipscHitz條件，既滿足存在常數(shù)H使 按照下式在和去倒數(shù)作線性組合（3）其中為待定系數(shù)，確定他們的準則是使（3）式有盡量

6、高的精度。注意到的假設，并對作二元泰勒展開，且利用，（3）式可寫為：于是故得截斷誤差為容易看出只要： 就可以使（3）式具有2階精度。 以上分析了龍格庫塔方法思想，下面用時的龍格庫塔方法即為改良后的歐拉公式法在MATHAB上解決本課程設計題目：3. .模型方程： 開始條件h=1.2m，d0=1.2m，g=9.8m*s(-2)，d=0.03m 水深為h時，流量Q為0.6(/4)d2*(gh)0.5，則水深下降dh所需時間 ：dt=-(/4)h2*dh/0.6(/4)d2*(gh)0.5=-h1.5*dh/0.6d2*(g)0.5水深由1.2m至0定積分得水從小孔流完的時間：T（其中已知d=0.03

7、m，g-9.8m/s） 設兩分鐘（120S）后水深為X m ，由 dt=-(/4)h2*dh/0.6*(/4)*d2*(gh)0.5=-h1.5*dh/0.6d2*(g)0.5 則263.93-120 =X2.5/1.5*d2*(g)0.5 以d=0.03，g=9.8代入上式得 水深：X 用歐拉方程和龍格庫塔方法則水深下降dh所需時間 ：dt=t(k+1)-t(k)=-(/4)h2*dh/0.6(/4)d2*(gh)0.5=-h1.5*dh/0.6d2*(g)0.5 k1=0.15*sqrt(g*(x(n)*d2/(-43.6359*x(n)8+213.0457*x(n)7-414.873*x

8、(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)*d2/(-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

9、)=x(n)-h*(k1+k2)/2;4. 編譯程序:(1)g=9.8;d=0.03;syms hy=-h(1.5)/(0.6*d2*sqrt(2*g);T=int(y,h,1.2,0);eval(T) x=(T-120)*(1.5*d2*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)*d2/(-43.6359*x(n)8+213.0457*x(n)7-

10、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)*d2/(-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,40