數(shù)值分析計(jì)算實(shí)習(xí)題

1、插值法下列數(shù)據(jù)點(diǎn)的插值x01491625364964y012345678可以得到平方根函數(shù)的近似，在區(qū)間0,64上作圖.用這9個(gè)點(diǎn)作8次多項(xiàng)式插值Ls(x).用三次樣條(第一邊界條件)程序求S(x).從得到結(jié)果看在0,64上，哪個(gè)插值更精確；在區(qū)間0,1上，兩種插值哪個(gè)更精確？解：(1)拉格朗日插值多項(xiàng)式，求解程序如下symsxl;x1=01491625364964;y1=012345678;n=length(x1);Ls=sym(0);fori=1:nl=sym(y1(i);fork=1:i-1l=l*(x-x1(k)/(x1(i)-x1(k);endfork=i+1:nl=l*(x-x1(

2、k)/(x1(i)-x1(k);endLs=Ls+l;endLs=simplify(Ls)%為所求插值多項(xiàng)式Ls(x).輸出結(jié)果為Ls=-24221063/63504000*x2+95549/72072*x-1/3048192000*x8-2168879/435456000*x4+19/283046400*x7+657859/10886400*x3+33983/152409600*x5T3003/2395008000*x飛(2)三次樣條插值，程序如下x1=01491625364964;y1=012345678;x2=0:1:64;y3=spline(x1,y1,x2);p=polyfit(x2

3、,y3,3);%得到三次樣條擬合函數(shù)S=p(l)+p(2)*x+p(3)*x2+p(4)*x3%得到S(x)輸出結(jié)果為：S=2288075067923491/73786976294838206464-2399112304472833/576460752303423488*x+4552380473376713/18014398509481984*x2+999337332656867/1125899906842624*x3在區(qū)間0,64上，分別對(duì)這兩種插值和標(biāo)準(zhǔn)函數(shù)作圖，plot(x2,sqrt(x2),b,x2,y2,r,x2,y3,y)藍(lán)色曲線為y=Vx函數(shù)曲線，紅色曲線為拉格朗日插值函數(shù)曲線

4、，黃色曲線為三次樣條插值曲線100806040200100806040200可以看到藍(lán)色曲線與黃色曲線幾乎重合，因此在區(qū)間0，64三次樣條插值更精確。在0，1區(qū)間上由上圖看不出差別，不妨代入幾組數(shù)據(jù)進(jìn)行比較，取x4二0:0.2:1x4=0:0.2:1;sqrt(x4)subs(Ls,x,x4)spline(x1,y1,x4)運(yùn)行結(jié)果為%準(zhǔn)確值%拉格朗日插值%三次樣條插值ans=00.44720.63250.77460.89441.0000ans=00.25040.47300.67060.84551.0000ans=00.24290.46300.66170.84031.0000從這幾組數(shù)值上可以

5、看出在0,1區(qū)間上，拉格朗日插值更精確。數(shù)據(jù)擬合和最佳平方逼近有實(shí)驗(yàn)給出數(shù)據(jù)表x0.00.10.20.30.50.81.0y1.00.410.500.610.912.022.46試求3次、4次多項(xiàng)式的曲線擬合，再根據(jù)數(shù)據(jù)曲線形狀，求一個(gè)另外函數(shù)的擬合曲線，用圖示數(shù)據(jù)曲線及相應(yīng)的三種擬合曲線。解：(1)三次擬合，程序如下symx;x0=0.00.10.20.30.50.81.0;y0=1.00.410.500.610.912.022.46;cc=polyfit(x0,y0,3);S3二cc(l)+cc(2)*x+cc(3)*x2+cc(4)*x3%三次擬合多項(xiàng)式xx=x0(1):0.1:x0(l

6、ength(x0);yy=polyval(cc,xx);plot(xx,yy,-);holdon;plot(x0,y0,x);xlabel(x);ylabel(y);運(yùn)行結(jié)果S3=-7455778416425083/1125899906842624+1803512222945437/140737488355328*x-655705280524945/140737488355328*x2+4172976892199509/4503599627370496*x3圖像如下(2)4次多項(xiàng)式擬合symx;x0=0.00.10.20.30.50.81.0;y0=1.00.410.500.610.912.0

7、22.46;cc=polyfit(x0,y0,4);S3二cc(l)+cc(2)*x+cc(3)*x2+cc(4)*x3+cc(5)*x4xx=x0(1):0.1:x0(length(x0);yy=polyval(cc,xx);plot(xx,yy,r);holdon;plot(x0,y0,x);xlabel(x);ylabel(y);運(yùn)行結(jié)果S3=3248542900396215/1125899906842624-3471944732519151/281474976710656*x+4580931990070637/281474976710656*x2-5965836931688425/n2

8、5899906842624*x3+8491275464650307/9007199254740992*x4圖像如下2.521.5y2.521.5y10.5x(3)另一個(gè)擬合曲線,新建一個(gè)M-file程序如下：functionC,L=lagran(x,y)w=length(x);n=w-1;L=zeros(w,w);fork=1:n+1V=1;forj=1:n+1ifk=jV=conv(V,poly(x(j)/(x(k)-x(j);endendL(k,:)=V;endC=y*L在命令窗口中輸入以下的命令：x=0.00.10.20.30.50.81.0;y=1.00.410.500.610.912

9、.022.46;cc=polyfit(x,y,4);xx=x(1):0.1:x(length(x);yy=polyval(cc,xx);plot(xx,yy,r);holdon;plot(x,y,x);xlabel(x);ylabel(y);x=0.00.10.20.30.50.81.0;y=0.940.580.470.521.002.002.46;%y中的值是根據(jù)上面兩種擬合曲線而得到的估計(jì)數(shù)據(jù)，不是真實(shí)數(shù)據(jù)C,L=lagran(x,y);xx=0:0.01:1.0;yy=polyval(C,xx);holdon;plot(xx,yy,b,x,y,.);圖像如下x解線性方程組的直接解法線性方

10、程組Ax=b的A及b為1078732A=75651078732A=7565，b=2386109337591031則解1111.用MATLAB內(nèi)部函數(shù)求detA及A的所有特征值和cond(A)2若令1078.17.2A+6A=708q504()65，求解(A+SA)(x+6x)=b，輸出向量x和85.989.8996.99599.98|6x|2從理論結(jié)果和實(shí)際計(jì)算兩方面分析線性方程組Ax=b解得相對(duì)誤差l|6x|2/|x|2及a的相對(duì)誤差|6a|2/|a|2的關(guān)系.解：(1)程序如下clear;A=10787;7565;86109;75910;det(A)cond(A，2)eig(A)輸出結(jié)果為

11、ans=1ans=2.9841e+003ans=0.01020.84313.858130.2887(2)程序如下A=1078.17.2;7.085.0465;85.989.899;6.99599.98;b=32233331;x0=1111;x=Ab%擾動(dòng)后方程組的解x1=x-x0norm(x1,2)x1=x-x0norm(x1,2)%8x的值%5x的2-范數(shù)運(yùn)行結(jié)果為x=-9.586318.3741-3.22583.5240 x1=-10.586317.3741-4.22582.5240ans=20.9322程序如下A=1078.17.2;7.085.0465;85.989.899;6.995

12、99.98;A0=10787;7565;86109;75910;b=32233331;x0=1111;x=Ab;x1=x-x0;norm(x1,2);A1=A-AO;%SA的值norm(xl,2)/norm(x0,2)%|x|2/|x|2的值norm(Al,2)/norm(A0,2)%|A|2/|A|2的值輸出結(jié)果為ans=10.4661ans=0.0076可見A相對(duì)誤差只為0.0076，而得到的結(jié)果x的相對(duì)誤差就達(dá)到了10.4661,該方程組是病態(tài)的，A的條件數(shù)為2984.1遠(yuǎn)遠(yuǎn)大于1,當(dāng)A只有很小的誤差就會(huì)給結(jié)果帶來很大的影響。非線性方程數(shù)值解法求下列方程的實(shí)根x八2-3x+2-e八x=0

13、;x八3+2x八2+10 x-20=0.要求：(1)設(shè)計(jì)一種不動(dòng)點(diǎn)迭代法，要使迭代序列收斂，然后再用斯特芬森加速迭代，計(jì)算到|x(k)-x(k-1)|10,-8)為止。(2)用牛頓迭代，同樣計(jì)算到|x(k)-x(k-1)|10,-8)。輸出迭代初值及各次迭代值和迭代次數(shù)k,比較方法的優(yōu)劣。解:(1)先用畫圖的方法估計(jì)根的范圍ezplot(，x“2_3*x+2_exp(x)，);gridon;可以估計(jì)到方程的根在區(qū)間(0,1)；選取迭代初值為x0=0.5；構(gòu)造不動(dòng)點(diǎn)迭代公式x(k+1)=(x(k廠2+2-e八x(k)/3;W(x)二(x八2+2-e八x)/3;當(dāng)0 x1時(shí)，0W(x)1;曠(x)

14、、1;因此迭代公式收斂。0.257534913615250.25753491361525程序如下：formatlong;f二inline(x2+2-exp(x)/3)disp(x=);x=feval(f,0.5);disp(x);Eps=1E-8;i=1;while1x0=x;i=i+1;x=feval(f,x);disp(x);ifabs(x-x0)Epsbreak;endendi,x運(yùn)行結(jié)果為f=Inlinefunction:f(x)=(x2+2-exp(x)/3x=0.200426243099960.272749065098370.253607156584130.258550376264

15、940.257265636335090.257598985162190.257512454514830.257529084167960.257530285439860.257530285439860.257530597238330.257530204510460.257530306445640.257530279987670.25753028685501i=14x=0.25753028685501斯特芬森加速法，程序如下:formatlong;f二inline(x-(x2+2-exp(x)/3-x)2/(x2+2-exp(x)/3)2+2-exp(x2+2-exp(x)/3)/3-2*(x2+

16、2-exp(x)/3+x);disp(x=);x=feval(f,0.5);disp(x);Eps=1E-8;i=1;while1x0=x;i=i+1;x=feval(f,x);disp(x);ifabs(x-x0)Epsbreak;endendi,x運(yùn)行結(jié)果為x=0.258684427565790.257530317719810.25753028543986i=4i=40.25753028543986i=40.25753028543986i=4x=0.25753028543986牛頓迭代法，程序如下：formatlong;x=sym(x);f二sym(，x“2-3*x+2-exp(x);df

17、=diff(f,x);FX=x-f/df;Fx=inline(FX);disp(，x=，);x1=0.5;disp(x1);Eps=1E-8;i=0;while1x0=x1;i=i+1;x1=feval(Fx,x1);disp(x1);ifabs(x1-x0)1E10break;endifabs(x-x0)Epsbreak;endendi,x運(yùn)行結(jié)果：f=Inlinefunction:f(x)=(-2*x2-10*x+20廠1/3x=0.166666666666676.09259259259259-38.38843164151806-8.478196837919431e+002-4.76366

18、0785374071e+005-1.512815059604763e+011i=6x=-1.512815059604763e+011迭代6次后x的值大得令人吃驚,表明構(gòu)造的式子并不收斂.也無法構(gòu)造出收斂的不動(dòng)點(diǎn)公式牛頓迭代法，程序如下：formatlong;x=sym(x);f=sym(x3+2*x2+10*x-20);df=diff(f,x);FX=x-f/df;Fx=inline(FX);disp(x=);x1=0.5;disp(x1);Eps=1E-8;i=0;while1x0=x1;i=i+1;x1=feval(Fx,x1);disp(x1);ifabs(x1-x0)Epsbreak;

19、endendi,x1運(yùn)行結(jié)果：x=1.500000000000001.373626373626371.368814819623961.368808107834411.36880810782137i=4x1=1.36880810782137比較三種方法，牛頓法的收斂性比較好，相比不動(dòng)點(diǎn)迭代法要構(gòu)造出收斂的公式比較難，牛頓法迭代次數(shù)也較少，收斂速度快，只是對(duì)初值的要求很高，幾種方法各有利弊，具體采用哪種也需因題而異。常微分方程初值問題數(shù)值解法給定初值問題y=-50y+50 x八2+2x,0 x1;y(0)=1/3;用經(jīng)典的四階R-K方法解該問題，步長分別取h=0.1,0.025,0.01,計(jì)算并打

20、印x=0.1i(i=0,l,，10)各點(diǎn)的值，與準(zhǔn)確值y(x)=l/3e八(-50 x)+x八2比較。解：取步長h=0.1,程序如下：%經(jīng)典的四階R-K方法clear;F二-50*y+50*x2+2*x;a=0;b=1;h=0.1;n=(b-a)/0.1;X=a:0.1:b;Y=zeros(1,n+1);Y(1)=1/3;fori=1:nx=X(i);y=Y(i);K1=h*eval(F);x=x+h/2;y=y+K1/2;K2=h*eval(F);y=Y(i)+K2/2;K3=h*eval(F);x=X(i)+h;y=Y(i)+K3;K4=h*eval(F);Y(i+1)=Y(i)+(K1+

21、2*K2+2*K3+K4)/6;end%準(zhǔn)確值temp=;f二dsolve(Dy=-50*y+50*x2+2*x,y(0)=l/3,x);df=zeros(1,n+1);fori=1:n+1temp=subs(f,x,X(i);df(i)=double(vpa(temp);enddisp（步長四階經(jīng)典R-K法準(zhǔn)確值）disp（X,Y,df）;運(yùn)行結(jié)果：步長四階經(jīng)典R-K法準(zhǔn)確值1.0e+010*00.000000000010000.0000000000200000.000000000010000.000000000020000.000000000030000.000000000040000.0

22、00000000050000.000000000060000.000000000070000.000000000080000.000000000090000.00000000010000%畫圖觀察結(jié)果figure;plot(X,df,k*,X,Y,gridon;0.000000000033330.000000000460550.000000006306250.000000086404940.000001184363000.000016235451100.000222560671340.003050935427780.041823239217400.573326903478097.8593563

23、00837710.000000000033330.000000000001220.000000000004000.000000000009000.000000000016000.000000000025000.000000000036000.000000000049000.000000000064000.000000000081000.00000000010000titletitle（四階經(jīng)典R-K法解常微分方程）;legend（準(zhǔn)確值，四階經(jīng)典R-K法）;當(dāng)X值接近1的時(shí)候，偏離準(zhǔn)確值太大。當(dāng)步長h=0.025時(shí)，將上面程序中的h改為0.025即可，運(yùn)行結(jié)果:步長00.10000000000

24、000步長00.100000000000000.200000000000000.300000000000000.400000000000000.500000000000000.600000000000000.700000000000000.800000000000000.900000000000001.00000000000000圖像如下：四階經(jīng)典R-K法0.333333333333330.103135240342880.044285273275990.051967957555070.093957311494390.160345314357570.248085701305570.356241884727580.484525906616270.6