數(shù)值分析上機(jī)實(shí)習(xí)報(bào)告

（數(shù)值分析上機(jī)實(shí)驗(yàn)報(bào)告）院系：礦業(yè)學(xué)院專(zhuān)業(yè)：礦業(yè)工程班級(jí)：2015姓名：王學(xué)號(hào)：2015022指導(dǎo)教師：代第一題1.用Newton法求解方程，在（0.1，1.9）的近似根（初始近似值取為區(qū)間端點(diǎn)，迭代6次或誤差小于0.00001）。1.1理論依據(jù)及方法應(yīng)用條件Newton迭代法:由一般迭代函數(shù)，取s=2時(shí)，有，可得二階迭代序列，此種迭代法稱(chēng)為Newton迭代法。定理:設(shè)函數(shù)在有限區(qū)間[a,b]上二階導(dǎo)數(shù)存在，且滿(mǎn)足條件（Ⅱ）在區(qū)間[a,b]上不變號(hào);（Ⅲ）≠0;（Ⅳ）||/b-a≤||其中c是a,b中使min[|,]達(dá)到的一個(gè);則對(duì)任意時(shí)近似值x0∈[a,b],由Newton迭代過(guò)程有:k=0,1,2…所產(chǎn)生的迭代序列{x0}平方收斂于方程=0區(qū)間[a,b]上的唯一解α。推論:設(shè)函數(shù)f(x)滿(mǎn)足定理中條件Ⅰ,Ⅱ,Ⅲ,若選初值，使·>0,則Newton迭代過(guò)程(k=0,1,2…)產(chǎn)生的迭代序列{xk}單調(diào)收斂于=0的唯一解α。1.2計(jì)算程序#include<iostream.h>#include<iomanip.h>#include<string>#include<math.h>usingnamespacestd;double*newton(doublea,doubleb,doubleeps); //牛頓迭代函數(shù)doublenewtonz(doublex); //牛頓迭代子函數(shù)voidmain() doublea=0.1,b=1.9,eps=0.00001,*result; //初始數(shù)據(jù) cout<<"\n牛頓法解方程：x^7-28x^4+14=0,在(0.1,1.9)中求近似根，初始值為區(qū)間端點(diǎn)，\n誤差為0.00001。\n"<<endl; cout<<"學(xué)號(hào)：2014021966姓名：徐林\n"<<endl; result=newton(a,b,eps); if(a<=result[0]&&result[0]<=b) cout<<"近似根為："<<result[0]<<endl; if(a<=result[1]&&result[1]<=b) cout<<"近似根為："<<result[1]<<endl; cout<<"\n"<<"結(jié)束,按任意鍵關(guān)閉"<<endl; getchar();}//主函數(shù)結(jié)束doublenewtonz(doublex) //牛頓迭代子函數(shù) doublex1=0.0,t; t=(7*pow(x,6)-4*28*pow(x,3)); if(t==0) exit(0); x1=(x-((pow(x,7)-28*pow(x,4)+14)/t)); returnx1;double*newton(doublea,doubleb,doubleeps) //牛頓迭代函數(shù) doublex0=0.0,x1=1.0,x2=0.0,re[2]; intk=0; x0=a; while(x0>eps) //代入a迭代計(jì)算 k++; x2=x1; x1=newtonz(x1); //調(diào)用牛頓迭代子函數(shù) x0=fabs(x1-x2); }re[0]=x1; x0=b,k=0; while(x0>eps) //代入b迭代計(jì)算 k++; x2=x1; x1=newtonz(x1); //調(diào)用牛頓迭代子函數(shù) x0=fabs(x1-x2); }re[1]=x1; returnre;1.3計(jì)算結(jié)果打印1.4MATLAB上機(jī)程序functiony=Newton(f,df,x0,eps,M)d=0;fork=1:Miffeval(df,x0)==0d=2;breakelsex1=x0-feval(f,x0)/feval(df,x0);ende=abs(x1-x0);x0=x1;ife<=eps&&abs(feval(f,x1))<=epsd=1;breakendendifd==1y=x1;elseifd==0y='迭代M次失敗';elsey='奇異'endfunctiony=df(x)y=7*x^6-28*4*x^3;Endfunctiony=f(x)y=x^7-28*x^4+14;End>>x0=1.9;>>eps=0.00001;>>M=100;>>x=Newton('f','df',x0,eps,M);>>vpa(x,7)1.5問(wèn)題討論1.需注意的是，要使用Newton迭代法須滿(mǎn)足定理中的條件Ⅰ,Ⅱ,Ⅲ，以及·>0。要用誤差范圍來(lái)控制循環(huán)的次數(shù)，保證循環(huán)的次數(shù)和質(zhì)量，編寫(xiě)程序過(guò)程中要注意標(biāo)點(diǎn)符號(hào)的使用，正確運(yùn)用適當(dāng)?shù)臉?biāo)點(diǎn)符號(hào)，Newton迭代法是局部收斂的，在使用時(shí)應(yīng)先確定初始值，否則所得的解可能不在所要求的范圍內(nèi)。(3)因?yàn)閚ewton法求方程是平方收斂的，所以較為精確,但是要求出函數(shù)的導(dǎo)數(shù)，且必須有二階導(dǎo)數(shù)。

第二題2.已知函數(shù)值如下表:1234500.693147181.09861231.38629441.60943786789101.79175951.94591012.0794452.19722462.3025851=1=0.1試用三次樣條插值求及的近似值。2.1理論依據(jù)及方法應(yīng)用條件三次樣條插值函數(shù)可定義為：對(duì)于[a,b]上的一個(gè)劃分∏a<x0<x1<x2<…<xn-1<xn=b.(n>=2)如果定義在[a,b]上的函數(shù)S(x),滿(mǎn)足（1）.在[xi,xi+1]上為3次多項(xiàng)式；（2）.S(x),S＇(x)，S＂(x)在[a,b]上連續(xù)，則稱(chēng)S(x)在[a,b]上劃分的3次樣條函數(shù)，如果對(duì)于，還滿(mǎn)足，，則稱(chēng)為的三次樣條插值函數(shù)。其基本思想是對(duì)均勻分劃的插值函數(shù)的構(gòu)造,三次樣條函數(shù)空間中1,x,,x2,x3,(x-xj)+3為基函數(shù),而取B樣條函數(shù)Ω3(（x-xj）/h)為基函數(shù).由于三次樣條函數(shù)空間是N+3維的,故我們把分點(diǎn)擴(kuò)大到X-1,XN+1,則任意三次樣條函數(shù)可用Ω3(（x-xj）/h)線性組合來(lái)表示S(x)=cjΩ3(（x-xj)/h)這樣對(duì)不同插值問(wèn)題,若能確定cj由解的唯一性就能求得S(x)。由s(xi)=yi,I=1,2,…Ns’(x0)=y0’,s’(xN)=yN’可得S(xi)=cjΩ3(（xi-xj）/h)=yiS’(x0)=1/hcjΩ3’(（x0-xj）/h)=y’0S’(xN)=1/hcjΩ3’(（xN-xj）/h)=y’N2.2計(jì)算程序#include<stdio.h>#include<math.h>#defineN10/*宏定義*/main()floats,ds,t;floatdy0=1,dy9=0.1;intj;intx[N]={1,2,3,4,5,6,7,8,9,10};floaty[N]={0,0.69314718,1.0986123,1.3862944,1.6094378,1.7917595,1.9459101,2.079445,2.1972246,2.3025851};intb[N]={2,2,2,2,2,2,2,2,2,2},h[N-1];floatd[N],u[N-1],v[N-1],a[N-1],c[N-1],B[N],l[N-1],p[N],X[N];for(j=1;j<=9;j++)h[j-1]=x[j]-x[j-1];d[0]=6/h[0]*(y[1]/h[0]-y[0]/h[0]-dy0);d[9]=6/h[8]*(dy9-y[9]/h[8]+y[8]/h[8]);for(j=1;j<=8;j++)d[j]=6/(h[j-1]+h[j])*(y[j+1]/h[j]-y[j]/h[j]-y[j]/h[j-1]+y[j-1]/h[j-1]);for(u[8]=1,j=0;j<=7;j++)u[j]=h[j-1]/(h[j-1]+h[j]);for(v[0]=1,j=1;j<=8;j++)v[j]=h[j]/(h[j-1]+h[j]);for(j=0;j<=8;j++)a[j]=u[j];for(j=0;j<=8;j++)c[j]=v[j];for(B[0]=b[0],j=1;j<=9;j++)/*追趕法求解三彎矩方程*/B[j]=b[j]-a[j]/B[j-1]*c[j-1];for(j=1;j<=9;j++)l[j]=a[j]/B[j-1];for(j=1;j<=9;j++)p[j]=d[j]-l[j]*p[j-1];X[9]=p[9]/B[9];for(j=8;j>=0;j--)X[j]=p[j]/B[j]-c[j]*X[j+1]/B[j];t=4.563;s=X[3]*pow((x[4]-t),3)/6/h[3]+X[4]*pow((t-x[3]),3)/6/h[3]+(y[3]-X[3]*h[3]*h[3]/6)*(x[4]/h[3]-t/h[3])+(y[4]-X[4]*h[3]*h[3]/6)*(t/h[3]-x[3]/h[3]);/*解f(x)的值*/ds=-X[3]*pow((x[4]-t),2)/2/h[3]+X[4]*pow((t-x[3]),2)/2/h[3]-(y[3]-X[3]*h[3]*h[3]/6)/h[3]+(y[4]-X[4]*h[3]*h[3]/6)/h[3];/*解f’(x)的值*/printf("s=%f\nds=%f\n",s,ds);/*打印結(jié)果*/2.3計(jì)算結(jié)果打印2.4MATLAB上機(jī)程序functionQ=san(ssss,p)Q=zeros(2,1);x=[1;2;3;4;5;6;7;8;9;10];y=[0;0.69314718;1.0986123;1.3862944;1.6094378;1.7917595;1.9459101;2.079445;2.1972246;2.3025851];h=zeros(10,1);d=zeros(10,1);u=zeros(10,1);v=zeros(10,1);r=zeros(10,1);l=zeros(10,1);z=zeros(10,1);m=zeros(10,1);fort=1:1:9;h(t)=x(t+1)-x(t);endd(1)=6/h(1)*((y(2)-y(1))/h(1)-1);d(10)=6/h(9)*(0.1-(y(10)-y(9))/h(9));fort=1:1:8u(t+1)=h(t)/(h(t)+h(t+1));v(t+1)=1-u(t+1);d(t+1)=6/(h(t)+h(t+1))*((y(t+2)-y(t+1))/(x(t+2)-x(t+1))-(y(t+1)-y(t))/(x(t+1)-x(t)));endu(10)=1;v(1)=1;r(1)=d(1);fort=2:1:10l(t)=u(t)/r(t-1);r(t)=d(t)-l(t)*v(t-1);endz(1)=d(1);fort=2:1:10z(t)=d(t)-l(t)*z(t-1);endm(10)=z(10)/r(10);fort=9:-1:1m(t)=(z(t)-v(t)*m(t+1))/r(t);endfort=1:1:10ifp>=t&&p<(t+1)Q(:,1)=feval(ssss,p,t,x,m,h,y);breakendendfunctionQ=ssss(p,t,x,m,h,y)Q=zeros(2,1);Q(1,1)=((power((x(t+1)-p),3)*m(t)+power((p-x(t)),3)*m(t+1))/6+(y(t)-m(t)*h(t)*h(t)/6)*(x(t+1)-p)+(y(t+1)-m(t+1)*h(t)*h(t)/6)*(p-x(t)))/h(t);Q(2,1)=(-(power((x(t+1)-p),2)*m(t)+power((p-x(t)),2)*m(t+1))/2+(y(t)-m(t)*h(t)*h(t)/6)+(y(t+1)-m(t+1)*h(t)*h(t)/6))/h(t);end2.5問(wèn)題討論1.若要用追趕法求解三對(duì)角方程組，三對(duì)角陣需要滿(mǎn)足：(i=1,2,…,n)均非奇異，保證有唯一的Doolittle分解;≠0;2.樣條插值效果比Lagrange插值好,三次樣條插值的解存在且唯一，近似誤差較小.并且沒(méi)有Runge現(xiàn)象。

第三題3.用Romberg算法求（允許誤差ε=0.00001）。3.1理論依據(jù)及方法應(yīng)用條件數(shù)值積分的Romberg算法計(jì)算步驟如下:當(dāng)時(shí)，就停機(jī)。3.2計(jì)算程序#include<stdio.h>#include<math.h>#defineN9floatf(floatx)/*定義函數(shù)f(x)*/floaty;y=pow(3,x)*pow(x,1.4)*(5*x+7)*sin(x*x);return(y);main()floatT[N+1][N+1],h[N+1],a=1,b=3,m[N+1];inti,l;T[1][0]=(b-a)*(f(a)+f(b))/2;l=1;while(l<=N)m[l]=0;for(i=1;i<=(pow(2,l-1));i++)m[l]+=f(a+(2*i-1)*(b-a)/pow(2,l));T[1][l]=(T[1][l-1]+(b-a)*m[l]/pow(2,l-1))/2;l++;i=1;while(i<=N)for(l=1;l<=N-i+1;l++)T[i+1][l-1]=(pow(4,i)*T[i][l]-T[i][l-1])/(pow(4,i)-1);h[i]=T[i][0]-T[i+1][0];if(fabs(h[i])<=1e-5)break;i++;printf("Theansweris:%f",T[i+1][0]);3.3計(jì)算結(jié)果打印3.4MATLAB上機(jī)程序function[T,n]=mromb(f,a,b,eps)ifnargin<4,eps=1e-6;endh=b-a;R(1,1)=(h/2)*(feval(f,a)+feval(f,b));n=1;J=0;err=1;while(err>eps)J=J+1;h=h/2;S=0;fori=1:nx=a+h*(2*i-1);S=S+feval(f,x);endR(J+1,1)=R(J,1)/2+h*S;fork=1:JR(J+1,k+1)=(4^k*R(J+1,k)-R(J,k))/(4^k-1);enderr=abs(R(J+1,J+1)-R(J+1,J));n=2*n;endR;T=R(J+1,J+1);formatlongf=@(x)(3.^x)*(x.^1.4)*(5*x+7)*sin(x*x);[T,n]=mromb(f,1,3,1.e-5)3.5問(wèn)題討論1、Romberge算法的優(yōu)點(diǎn)是:a、把積分化為代數(shù)運(yùn)算,而實(shí)際上只需求T1(i),以后用遞推可得。b、算法簡(jiǎn)單且收斂速度快,一般4或5次即能達(dá)到要求。c、節(jié)省存儲(chǔ)量，算出的EMBEDEquation.3可存入EMBEDEquation.3。2、Romberge算法的缺點(diǎn)是:a、對(duì)函數(shù)的光滑性要求較高。b、計(jì)算新分點(diǎn)的值時(shí),這些數(shù)值的個(gè)數(shù)成倍增加。

第四題4.用定步長(zhǎng)四階Runge-Kutta法求解,打印,,,,4.1理論依據(jù)及方法應(yīng)用條件Runge-Kutta法的基本思想:不是按Taylor公式展開(kāi)，而是先寫(xiě)成處附近的值的線性組合（有待定系數(shù)）,再按Taylor公式展開(kāi)，然后確定待定常數(shù)。四階古典Runge-Kutta公式:4.2計(jì)算程序#include<stdio.h>intmain() inti; doubleh=0.0005; doublek1,k2,k3,k4; doubley1=0.0,y2=0.0,y3=0.0; for(i=1;i<=200;i++) k1=k2=k3=k4=h*1.0; y1+=(k1+2*k2+2*k3+k4)/6; k1=k2=k3=k4=h*y3; y2+=(k1+2*k2+2*k3+k4)/6; k1=h*(1000-1000*y2-100*y3); k2=h*(1000-1000*y2-100*(y3+0.5*k1)); k3=h*(1000-1000*y2-100*(y3+0.5*k2)); k4=h*(1000-1000*y2-100*(y3+k3)); y3+=(k1+2*k2+2*k3+k4)/6; if(i==50) printf("\ny1(0.025)=%fy2(0.025)=%fy3(0.025)=%f",y1,y2,y3); continue; if(i==90) printf("\ny1(0.045)=%fy2(0.045)=%fy3(0.045)=%f",y1,y2,y3); continue; if(i==170) printf("\ny1(0.085)=%fy2(0.085)=%fy3(0.085)=%f",y1,y2,y3); continue; if(i==200) printf("\ny1(0.100)=%fy2(0.100)=%fy3(0.100)=%f\n\n",y1,y2,y3);4.3計(jì)算結(jié)果打印4.4MATLAB上機(jī)程序functionY=R_K(df1,a,b,h)m=(b-a)/h;Y=zeros(3,1);S=zeros(3,1);K=zeros(3,4);x=a;y1=a;y2=a;y3=a;forn=1:mK(:,1)=feval(df1,x,y1,y2,y3);x=x+0.5*h;S(:,1)=Y+0.5*h.*K(:,1);y1=S(1,1);y2=S(2,1);y3=S(3,1);K(:,2)=feval(df1,x,y1,y2,y3);S(:,1)=Y+0.5*h.*K(:,2);y1=S(1,1);y2=S(2,1);y3=S(3,1);K(:,3)=feval(df1,x,y1,y2,y3);x=x+0.5*h;S(:,1)=Y+h.*K(:,3);y1=S(1,1);y2=S(2,1);y3=S(3,1);K(:,4)=feval(df1,x,y1,y2,y3);Y=Y+h.*(K(:,1)+2.*K(:,2)+2.*K(:,3)+K(:,4))/6;endfunctionZ=df1(x,y1,y2,y3)Z=zeros(3,1);Z(1)=0*x+0*y1+0*y2+0*y3+1;Z(2)=0*x+0*y1+0*y2+1*y3;Z(3)=0*x+0*y1-1000*y2-100*y3+1000;end4.5問(wèn)題討論1.定步長(zhǎng)四階runge-kutta法穩(wěn)定，精度高，可根據(jù)有變化的情況與需要的精度自動(dòng)修改步長(zhǎng)，誤差小且程序簡(jiǎn)單，存儲(chǔ)量少。2.但是Runge-Kutta法需要每步都計(jì)算函數(shù)值四次,在函數(shù)較復(fù)雜時(shí),工作量就會(huì)變得較大,可靠性有待核查。

第五題5.已知A與bA=12.384122.115237-1.0610741.112336-0.1135840.7187191.7423823.067813-2.0317432.11523719.141823-3.125432-1.0123452.1897361.563849-0.7841561.1123483.123124-1.061074-3.12543215.5679143.1238482.0314541.836742-1.0567810.336993-1.0101031.112336-1.0123453.12384827.1084374.101011-3.7418562.101023-0.71828-0.037585-0.1135842.1897362.0314544.10101119.8979180.431637-3.1112232.1213141.7843170.7187191.5638491.836742-3.7418560.4316379.789365-0.103458-1.1034560.2384171.742382-0.784165-1.0567812.101023-3.111223-0.10345814.7138463.123789-2.2134743.0678131.1123480.336993-0.718282.121314-1.1034563.12378930.7193344.446782-2.0317433.123124-1.010103-0.0375851.7843170.238417-2.2134744.44678240.00001b=(2.187436933.992318-25.1734170.846716951.784317-86.6123431.11012304.719345-5.6784392)T(1).用列主元素消去法求解Bx=b.5.1（1）理論依據(jù)用Househloder變換，把A化為三對(duì)角陣。設(shè)A=（aij）,aij=aji,A=(a1,a2,…,an)，ai=(a1i,a2i,…,ani)T,取第一列a1中a11不動(dòng)，把（a21,a31,…,an1）T化為（±s1,0,0,…,0）T,這里,前面的正負(fù)表示鏡像可取為兩個(gè)方向相反的向量，這增加計(jì)算的靈活性。記d1=(a11,±s1,0,…,0)T,則d1=H1a1,這里為了增加計(jì)算的靈活性，避免同號(hào)數(shù)相減，選取符號(hào)使2a21s1>0,這樣只需要把上式右端第二項(xiàng)取為2sign(a21)a21s1即可，從而歸納起來(lái)算法步驟為：（2）程序主體#include<stdio.h>#include<math.h>voidmain() inti,j,m,r,sign; doubleA0[9][9],s,z,p,n,v,h,l,y[9],u[9],k,q[9],X[9],x[9]={0,0,0,0,0,0,0,0,0},B[9][9],g[9]; doubleA[9][9]= {12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743}, {2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124}, {-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103}, {1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.037585}, {-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784317}, {0.718719,1.563849,1.836742,-3.741865,0.431637,9.789365,-0.103458,-1.103456,0.238417}, {1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.713846,3.123789,-2.213474}, {3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782}, {-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00001} doubleb[9]={2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392}; for(i=0;i<9;i++) for(j=0;j<9;j++) A0[i][j]=A[i][j]; for(r=0;r<7;r++) s=0; for(i=r+1;i<9;i++) s=s+A[i][r]*A[i][r]; s=sqrt(s); l=s*s+fabs(A[r+1][r])*s; if(A[r+1][r]>0)sign=1; elseif(A[r+1][r]<0)sign=-1; for(i=0;i<9;i++) if(i<=r)u[i]=0; elseif(i==r+1)u[i]=A[r+1][r]+sign*s; elseu[i]=A[i][r]; for(i=0;i<9;i++) y[i]=0; for(j=0;j<9;j++) y[i]=y[i]+A[i][j]*u[j]; y[i]=y[i]/l; k=0; for(i=0;i<9;i++) k=k+u[i]*y[i]; k=k/(2*l); for(i=0;i<9;i++) q[i]=y[i]-k*u[i]; for(i=0;i<9;i++) for(j=0;j<9;j++) A[i][j]=A[i][j]-q[i]*u[j]-q[j]*u[i];printf("\nHousehold變換矩陣B為:\n");for(i=0;i<9;i++) printf("\n"); for(j=0;j<9;j++) printf("%f,",A[i][j]);（3）結(jié)果截屏(4)MATLAB上機(jī)程序unction[x]=mgauss2(A,b,flag)ifnargin<3,flag=0;endn=length(b);fork=1:(n-1)[ap,p]=max(abs(A(k:n,k)));