南京工業(yè)大學(xué)工程數(shù)學(xué)第5章作業(yè)

1、學(xué)號(hào)：642081701046 姓名：苗潔 班級(jí)：化工2班第五章作業(yè)1.解：第一類(lèi)邊界條件的三次樣條插值計(jì)算如下：>>x=0.520,3.1,8.0,17.95,28.65,39.62,50.65,78,104.6,156.6,208.6,260.7,312.5,364.4,416.3,468,494,507,520;>>y=5.288,9.4,13.84,20.20,24.90,28.44,31.10,35,36.9,36.6,34.6,31.0,26.34,20.9,14.8,7.8,3.7,1.5,0.2;>> pp2=csape(x,y,'c

2、omplete',1.86548,-0.046115)pp2 = form: 'pp' breaks: 0.5200 3.1000 8 17.9500 28.6500 39.6200 50.6500 78 104.6000 156.6000 208.6000 260.7000 312.5000 364.4000 416.3000 468 494 507 520 coefs: 18x4 double pieces: 18 order: 4 dim: 1pp2.coefsans =-0.0029 -0.0979 1.8655 5.28800.0080 -0.1201 1.3

3、031 9.4000-0.0003 -0.0029 0.7007 13.84000.0003 -0.0128 0.5447 20.20000.0000 -0.0039 0.3655 24.90000.0000 -0.0038 0.2803 28.44000.0000 -0.0030 0.2052 31.1000-0.0000 -0.0009 0.0990 35.00000.0000 -0.0013 0.0403 36.9000-0.0000 -0.0000 -0.0301 36.60000.0000 -0.0004 -0.0527 34.60000.0000 -0.0002 -0.0816 3

4、1.00000.0000 -0.0001 -0.0979 26.3400-0.0000 -0.0001 -0.1109 20.9000-0.0000 -0.0002 -0.1254 14.8000-0.0000 -0.0002 -0.1462 7.80000.0001 -0.0009 -0.1748 3.70000.0000 0.0030 -0.1467 1.5000>> fnplt(pp2)再求x=2,30,130,350,515各點(diǎn)上的函數(shù)值，則輸入>> xi=2,30,130,350,515;>> yi=ppval(pp2,xi)yi =7.8252

5、25.3862 37.2138 22.4751 0.5427三次樣條插值函數(shù)為=+-+-由以上結(jié)果可得：當(dāng)0.520，3.1時(shí)：當(dāng)28.65，39.62時(shí)：當(dāng)104.6，156.6時(shí)：當(dāng)312.5，364.4時(shí)：當(dāng)507，520時(shí)：綜上，計(jì)算結(jié)果為：2301303505157.825225.386237.213822.47510.54274.11210.3296-0.1053-0.0991-0.05742.解：（1）復(fù)化公式的余項(xiàng)估計(jì),=35采用復(fù)化的要求具有9位有效數(shù)字則至少要取35個(gè)點(diǎn)，則利用MATLAB程序計(jì)算先求的四階導(dǎo)數(shù)：>> syms x f=4/(1+x2) %定義函

6、數(shù)f（x） n=input('輸入所求導(dǎo)數(shù)階數(shù):') f2=diff(f,x,n)f =4/(x2 + 1) 輸入所求導(dǎo)數(shù)階數(shù):4n = 4f2 =96/(x2+1)3-(1152*x2)/(x2+1)4+(1536*x4)/(x2+1)5復(fù)化辛普森程序：>> syms x f=inline('(4/(1+x2)','x');f2=inline('(96/(x2 + 1)3 - (1152*x2)/(x2 + 1)4 + (1536*x4)/(x2 + 1)5)','x'); f3='-(96/

7、(x2 + 1)3 - (1152*x2)/(x2 + 1)4 + (1536*x4)/(x2 + 1)5)' ; e=1*10(-9); a=0 ; b=1; x1=fminbnd(f3,0,1); for n=2:100000 Rn=-(b-a)/180*(b-a)/(2*n)4*f2(x1) ; if abs(Rn)<e break endendh=(b-a)/n ; Sn1=0 ; Sn2=0 ;for k=0:n-1 xk=a+k*h ; xk1=xk+h/2 ; Sn1=Sn1+f(xk1) ; Sn2=Sn2+f(xk) ;endSn=h/6*(f(a)+4*Sn1

8、+2*(Sn2-f(a)+f(b); z=exp(2); Rn=Sn-z ; format long;fprintf('Sn=') disp(Sn) fprintf('n=') disp(n) fprintf('Rn=') disp(Rn)計(jì)算結(jié)果：Sn= 3.141592653589789n= 76Rn= -4.247463445340861結(jié)果分析：與上述的余項(xiàng)估算結(jié)果對(duì)比可知，在其估算的范圍之內(nèi)，但是Rn的值不符合題目要求，不知錯(cuò)在哪里還請(qǐng)老師指正。（2）被積函數(shù)的文件M為：function y=fun553(x)y=4/(1+x2);Ro

9、mberg程序文件：function s,n=romberg(fun,a,b,n0,e)%Initial Values of Multiple-Trapezoidal integrationh=(b-a)/n0;sum=0;for i=1:n0-1 sum=sum+feval(fun,a+i*h);endr(1,1)=h*(feval(fun,a)+2*sum+feval(fun,b)/2;s1=r(1,1);s0=0;n=1;%Romberg Methodwhile(abs(s1-s0)>=e)|(h>=(b-a)/(24) sum=0; for j=1:2(n-1)*n0 su

10、m=sum+feval(fun,a+(j-0.5)*h); end n=n+1;h=h/2;r(n,1)=r(n-1,1)/2+sum*h; for j=2:n k=n-j+1; r(k,j)=r(k+1,j-1)+(r(k+1,j-1)-r(k,j-1)/(4(j-1)-1); end s1=r(1,n);s0=r(1,n-1);end%Compute resultsr;s=s1;n;調(diào)用該程序計(jì)算定積分：a=0;b=1;n0=1;e=1e-9;format long；s1,n=romberg('fun553',a,b,n0,e)，format short計(jì)算結(jié)果為：s1 =

11、3.141592653589723n =7（3）余項(xiàng)估計(jì)：=32Matlab程序：調(diào)用nagsint文件，然后輸入：>> format long;nagsint(inline('fun553'),eps,1,32,2),format short;ans = 0.245489451021956 0.254510548978044則其結(jié)果為0.5，但與其結(jié)果=3.1415926不相符。結(jié)果分析：這個(gè)計(jì)算結(jié)果是不正確的，我將老師給的程序和例子重新調(diào)用了一遍，其結(jié)果與老師給出的結(jié)果一致，說(shuō)明程序是正確的；在繼續(xù)找錯(cuò)誤的過(guò)程中，我將函數(shù)文件夾中的被積函數(shù)改成一個(gè)常數(shù)，其預(yù)算結(jié)

12、果也是正確的；我上網(wǎng)查閱了資料，也采用了很多方法，此函數(shù)最終的結(jié)果還是如此。3.解：（1）采用三次樣條插值計(jì)算：>>x1=0,0.921,1.843,2.949,3.871,4.978,5.900,7.006,7.928,8.967;x2=10.954,12.032,12.954,13.875,14.982,15.903,16.826,17.931,19.037,19.959,20.839;x3=22.958,23.880,24.986,25.908;>>y1=9.677,9.479,9.308,9.125,8.982,8.814,8.686,8.525,8.388,8

13、.220;y2=10.820,10.500,10.210,9.936,9.653,9.409,9.180,8.921,8.662,8.433,8.220;y3=10.820,10.591,10.354,10.180;>> pp1=spline(x1,y1)pp1 = form: 'pp' breaks: 0 0.9210 1.8430 2.9490 3.8710 4.9780 5.9000 7.0060 7.9280 8.9670 coefs: 9x4 double pieces: 9 order: 4 dim: 1>> pp2=spline(x2,y2

14、)pp2 = form: 'pp' breaks: 10.9540 12.0320 12.9540 13.8750 14.9820 15.9030 16.8260 17.9310 19.0370 19.9590 20.8390 coefs: 10x4 double pieces: 10 order: 4 dim: 1>> pp3=spline(x3,y3)pp3 = form: 'pp' breaks: 22.9580 23.8800 24.9860 25.9080 coefs: 3x4 double pieces: 3 order: 4 dim:

15、1>> pp1.coefs,pp2.coefs,pp3.coefsans = -0.0023 0.0224 -0.2336 9.6770 -0.0023 0.0160 -0.1983 9.4790 -0.0013 0.0097 -0.1746 9.3080 -0.0027 0.0055 -0.1578 9.1250 0.0042 -0.0020 -0.1547 8.9820 -0.0072 0.0120 -0.1437 8.8140 0.0029 -0.0081 -0.1401 8.6860 -0.0028 0.0014 -0.1475 8.5250 -0.0028 -0.0064

16、 -0.1521 8.3880ans = 0.0055 -0.0257 -0.2755 10.8200 0.0055 -0.0080 -0.3118 10.5000 0.0100 0.0071 -0.3126 10.2100 -0.0166 0.0349 -0.2739 9.9360 0.0136 -0.0203 -0.2578 9.6530 -0.0042 0.0173 -0.2605 9.4090 -0.0012 0.0058 -0.2392 9.1800 -0.0041 0.0017 -0.2310 8.9210 0.0056 -0.0118 -0.2423 8.6620 0.0056

17、0.0037 -0.2497 8.4330ans = -0.0014 0.0210 -0.2665 10.8200 -0.0014 0.0171 -0.2314 10.5910 -0.0014 0.0123 -0.1989 10.3540由這三組數(shù)據(jù)，可得區(qū)間8.388,8.967,10.954,12.032,22.958,23.880的插值多項(xiàng)式當(dāng)8.388,8.967時(shí)：當(dāng)10.954,12.032時(shí)：當(dāng)22.958,23.880時(shí)：因?yàn)檫@3組數(shù)是以觀測(cè)點(diǎn)為界的，可見(jiàn)這三點(diǎn)作為這三個(gè)函數(shù)的臨界點(diǎn)，通過(guò)求導(dǎo)，分別將9.9811，10.925，22.015帶入可得,(3) >> x

18、=9.9811,10.925,22.015;>> y=0.1929,0.2740,0.3098;>> pp1=spline(x,y)pp1 = form: 'pp' breaks: 9.9811 22.0150 coefs: -0.0069 0.0924 0.1929 pieces: 1 order: 3 dim: 1>> pp1.coefsans = -0.0069 0.0924 0.1929>> fnplt(pp1)可得水流速度的函數(shù)圖如下第（4）-（7）題，我按照老師所給的程序去執(zhí)行，總是出現(xiàn)無(wú)效表達(dá)“ppder”，通過(guò)上網(wǎng)

19、查找資料修改，結(jié)果還是錯(cuò)誤，還請(qǐng)老師幫我找出原因。第六章作業(yè)1.解：引入輔助變量，則可將原方程組化為一階方程組：，其中，輸入函數(shù)，function f=naeg(t,y)f=y(3);y(4);2*y(4)+y(1)-(a*y(1)+b)/c13-(b*(y(1)-a)/c23;-2*y(3)+y(2)-a*y(2)/c13-b*y(2)/c23;a=1-b,b=1/82.45;c1=(y(1)+b)2+y(2)2)0.5;c2=(y(1)+a)2+y(2)2)0.5;在命令窗口執(zhí)行>> y0=1.2,0,0,-1.04935371;>> t,y=ode45(naeg,

20、0,4,y0);t,yans = 0 1.2000 0 0 -1.0494 0.0000 1.2000 -0.0000 -0.0001 -1.0494 0.0001 1.2000 -0.0001 -0.0001 -1.0494 0.0001 1.2000 -0.0001 -0.0002 -1.0494 0.0001 1.2000 -0.0001 -0.0002 -1.0494 0.0003 1.2000 -0.0003 -0.0005 -1.0494 0.0004 1.2000 -0.0005 -0.0007 -1.0494 0.0006 1.2000 -0.0006 -0.0010 -1.0

21、494 0.0008 1.2000 -0.0008 -0.0012 -1.04940.0016 1.2000 -0.0016 -0.0025 -1.0494(數(shù)據(jù)太多，省略了中間的) 3.4461 3.2344 0.0306 15.1934 -2.9719 3.4916 3.9221 -0.1360 15.0195 -4.3495 3.5371 4.5998 -0.3651 14.7529 -5.7160 3.5826 5.2633 -0.6560 14.3937 -7.0658 3.6869 6.7090 -1.5512 13.2260 -10.0688 3.7913 8.0082 -2.7

22、509 11.5944 -12.8885 3.8956 9.1135 -4.2326 9.5252 -15.46104.0000 9.9817 -5.9674 7.0535 -17.7273>> plot(t,y(:,1),t,y(:,2),':')2.解：首先建立M函數(shù)function f=f_shoot(y,x)f=y(2);-y(1)*y(2)/4+x3/2+4;然后在命令窗口，輸入打靶法的程序>> e=0.5e-6; i=1; m(1)=0; m(2)=10; yn=35/3; h=0.02; N=1/h;while 1 x(1)=2; y(:,

23、1)=8;m(i); for n=1:N k1=f_shoot(y(:,n),x(n); k2=f_shoot(y(:,n)+h*k1/2,x(n)+h/2); k3=f_shoot(y(:,n)+h*k2/2,x(n)+h/2); k4=f_shoot(y(:,n)+h*k3,x(n)+h); y(:,n+1)=y(:,n)+h*(k1+2*k2+2*k3+k4)/6; x(n+1)=x(n)+h; end y1_end(i)=y(1,n+1); if (abs(y1_end(i)-yn)<e) x',y', break, end if i>=2 m(i+1)=m

24、(i)-(m(i)-m(i-1)*(y1_end(i)-yn)/(y1_end(i)-y1_end(i-1); end i=i+1;endans = 2.0000 8.0000 2.0000 2.0200 8.0408 2.0794 2.0400 8.0832 2.1577 2.0600 8.1271 2.2348 2.0800 8.1726 2.3109 2.1000 8.2195 2.3859 2.1200 8.2680 2.4600 2.1400 8.3179 2.5331 2.1600 8.3693 2.6053 2.1800 8.4221 2.6766 2.2000 8.4764 2

25、.7471 2.2200 8.5320 2.8168 2.2400 8.5890 2.8856 2.2600 8.6474 2.9537 2.2800 8.7072 3.0211 2.3000 8.7683 3.0877 2.3200 8.8307 3.1537 2.3400 8.8944 3.2190 2.3600 8.9594 3.2836 2.3800 9.0257 3.3477 2.4000 9.0933 3.4111 2.4200 9.1622 3.4740 2.4400 9.2323 3.5363 2.4600 9.3036 3.5980 2.4800 9.3762 3.6593

26、2.5000 9.4500 3.7200 2.5200 9.5250 3.7802 2.5400 9.6012 3.8400 2.5600 9.6786 3.8993 2.5800 9.7572 3.9581 2.6000 9.8369 4.0166 2.6200 9.9178 4.0746 2.6400 9.9999 4.1322 2.6600 10.0831 4.1894 2.6800 10.1675 4.2462 2.7000 10.2530 4.3026 2.7200 10.3396 4.3587 2.7400 10.4273 4.4144 2.7600 10.5162 4.4698

27、2.7800 10.6061 4.5249 2.8000 10.6971 4.5796 2.8200 10.7893 4.6340 2.8400 10.8825 4.6881 2.8600 10.9768 4.7420 2.8800 11.0722 4.7955 2.9000 11.1686 4.8488 2.9200 11.2661 4.9017 2.9400 11.3647 4.9545 2.9600 11.4643 5.0069 2.9800 11.5650 5.0591 3.0000 11.6667 5.1111然后繪圖：>> plot(x,y(1,:),'-r&#

28、39;,x,y(2,:),'.-b');xlabel('x'),ylabel('y(x) and yprime(x)');text(2.1,7,'Guess for yprime(0)=',num2str(m(i);text(2.6,7,'y(3)=',num2str(y1_end(i);gtext('y(x)');gtext('yprime(x)');m0=m，由題意可知藍(lán)線是，紅線為3.解：基于四階RK公式的打靶法求解桿沿軸x方向的溫度分布首先建立M函數(shù)function f=f_

29、shoot(y,x,a,b)f=y(2);a*(y(1)-b);然后在命令窗口，輸入打靶法的程序>> e=1e-6; i=1; m(1)=0; m(2)=100; r1=0;while 1 A=1e-4; p=1e-2; hc=20; k=60; b=293; a=p*hc/A/k; h=0.01; N=0.2/h; x(1)=0; y(:,1)=493;m(i); for n=1:N k1=f_shoot(y(:,n),x(n),a,b); k2=f_shoot(y(:,n)+h*k1/2,x(n)+h/2,a,b); k3=f_shoot(y(:,n)+h*k2/2,x(n)+h/2,a,b); k4=f_shoot(y(:,n)+h*k3,x(n)+h,a,b); y(:,n+1)=y(:,n)+h*(k1+2*k2+2*k3+k4)/6; x(n+1)=n*h; end y2_end(i)=y(2,n+1); if (abs(y2_end(i)-r1)<e) x