北航数值分析大作业题目三
北航数值分析大作业3
数值分析第三次作业1.设计方案对Fredholm积分方程，用迭代法进行求解：()'(())u x A u x=，其中11(())()(,)()A u x g x K x y u y dy-=-⋅⎰对于公式中的积分部分用数值积分方法。
复化梯形积分法，取2601个节点，取迭代次数上限为50次。
实际计算迭代次数为18次，最后算得误差为r= 0.97E-10。
复化Simpson积分法，取迭代次数上限为50次，取2*41+1，即83个节点时能满足精度要求。
实际计算迭代次数为17次，最后的误差为r= 0.97E-10。
Guass积分法选择的Gauss—Legendre法，取迭代次数上限为50次，直接选择8个节点，满足精度要求。
实际计算迭代次数为24次，最后算得误差为r= 0.87E-10。
2.全部源程序module integralimplicit nonecontains!//////////复化梯形subroutine trapezoid(m)implicit noneinteger :: i,j,k,mreal*8 :: x(m+1),u(m+1)real*8 :: sum,sum1,g,r,hreal*8 :: e=1.0e-10h=2./mdo i=1,m+1x(i)=-1.+(i-1)*hend dou=0.02do k=1,50do i=1,m+1sum1=0.g=dexp(x(i)*4.)+(dexp(x(i)+4.)-dexp(-4.-x(i)))/(x(i)+4.)do j=2,msum1=sum1+dexp(x(i)*x(j))*u(j)end dosum=h/2.*(dexp(x(i)*-1.)*u(1)+dexp(x(i)*1.)*u(m+1)+2*sum1)u(i)=g-sumend dor=h/2.*((dexp(x(1)*4)-u(1))**2+(dexp(x(m+1)*4)-u(m+1))**2) do i=2,mr=r+h*(dexp(x(i)*4)-u(i))**2end doif(dabs(r)<=e) exitend dowrite(*,*) kopen(1,file="trapezoid.txt")do i=1,m+1write(1,'(3(f18.12))') x(i),u(i),dexp(x(i)*4.)end dowrite(1,'(4x,a2,e9.2)') "r=",rclose(1)returnend subroutine trapezoid!///////////复化simpsonsubroutine simpson(m)implicit noneinteger :: i,j,k,mreal*8 :: x(2*m+1),u(2*m+1)real*8 :: sum,sum1,sum2,g,r,hreal*8 :: e=1.0e-10h=2./(2.*m)do i=1,2*m+1x(i)=-1.+(i-1)*hend dou=0.02do k=1,50do i=1,2*m+1sum1=0.sum2=0.g=dexp(x(i)*4.)+(dexp(x(i)+4.)-dexp(-4.-x(i)))/(x(i)+4.)do j=1,msum1=sum1+dexp(x(i)*x(2*j))*u(2*j)end dodo j=1,m-1sum2=sum2+dexp(x(i)*x(2*j+1))*u(2*j+1)sum=h/3.*(dexp(x(i)*-1.)*u(1)+dexp(x(i)*1.)*u(2*m+1)+4*sum1+2*sum2) u(i)=g-sumend dor=h/3.*((dexp(x(1)*4)-u(1))**2+(dexp(x(2*m+1)*4)-u(2*m+1))**2)do i=1,mr=r+4.*h/3.*(dexp(x(2*i)*4)-u(2*i))**2end dodo i=1,m-1r=r+2.*h/3.*(dexp(x(2*i+1)*4)-u(2*i+1))**2end doif(dabs(r)<=e) exitend dowrite(*,*) kopen(2,file="simpson.txt")do i=1,2*m+1write(2,'(3(f18.12))') x(i),u(i),dexp(x(i)*4.)end dowrite(2,'(4x,a2,e9.2)') "r=",rclose(2)returnend subroutine simpson!///////////Gauss_Legendre法subroutine Gaussimplicit noneinteger,parameter :: m=8integer :: i,j,kreal*8 :: x(m),u(m),a(m)real*8 :: sum,g,rreal*8 :: e=1.0e-10data x /-0.9602898565,-0.7966664774,-0.5255324099,-0.1834346425,&0.1834346425,0.5255324099,0.7966664774,0.9602898565/data a /0.1012285363,0.2223810345,0.3137066459,0.3626837834,&0.3626837834,0.3137066459,0.2223810345,0.1012285363/u=0.02do k=1,50do i=1,mg=dexp(x(i)*4.)+(dexp(x(i)+4.)-dexp(-4.-x(i)))/(x(i)+4.)do j=1,msum=sum+dexp(x(i)*x(j))*u(j)*a(j)end dou(i)=g-sumend dor=0.do i=1,mr=r+a(i)*(dexp(x(i)*4)-u(i))**2end doif(dabs(r)<=e) exitend dowrite(*,*) kopen(3,file="Gauss.txt")do i=1,mwrite(3,'(3(f18.12))') x(i),u(i),dexp(x(i)*4.)end dowrite(3,'(4x,a2,e9.2)') "r=",rclose(3)returnend subroutine Gaussend module!//////////主程序program mainuse integralimplicit noneinteger :: code1=2600integer :: code2=41call trapezoid(code1)call simpson(code2)call Gaussend program3.各种积分方法的节点和数值解（由于数据太多，在打印时用了较计算时少的有效数字）复化Simpson法4.各方法所得曲线（由于所取节点太多，且精度高，所以图中很难看出各曲线的区别。
北航研究生数值分析上机作业 三 (报告+所有程序大全)
数值分析上机作业3——求解非线性方程组以及二元函数的插值拟合1. 算法设计对于全部的插值节点(,),0,1,...,10,0,1,...,20i j x y i j ==，带入非线性方程组中，用Newton 迭代法解非线性方程组，得到(,),0,1,...,10,0,1,...,20i j t u i j ==。
对(,)i j t u ，在二维数表中进行插值，采用分片双二次插值法。
插值过程中，先选择分片区域的中心节点，在数表中的列记为(0:5)tt ，行记为(0:5)uu ，中心节点记为(,)a b ，生成向量_(0:2)t temp ，_(0)(())((1))/(((1)())((1)(1)))i i t temp t tt a t tt a tt a tt a tt a tt a =--+----+， _(1)((1))((1))/((()(1))(()(1)))i i t temp t tt a t tt a tt a tt a tt a tt a =---+---+， _(2)((1))(())/(((1)(1))((1)()))i i t temp t tt a t tt a tt a tt a tt a tt a =---+--+-，同理，生成向量_(0:2)u temp ，_(0)(())((1))/(((1)())((1)(1)))_(1)((1))((1))/((()(1))(()(1)))_(2)((1))(())/(((1)(1))((1)())j j j j j j u temp u uu a u uu a uu a uu a uu a uu a u temp u uu a u uu a uu a uu a uu a uu a u temp u uu a u uu a uu a uu a uu a uu a =--+----+=---+---+=---+--+-)记数表中以分片区域中心节点为中心的3×3的矩阵为T , 对于(,)i j t u 插值结果为(_)()(_)T t temp T u temp 。
数值分析大作业三、四、五、六、七
大作业 三1. 给定初值0x 及容许误差，编制牛顿法解方程f (x )=0的通用程序.解：Matlab 程序如下：函数m 文件：function Fu=fu(x) Fu=x^3/3-x; end函数m 文件：function Fu=dfu(x) Fu=x^2-1; end用Newton 法求根的通用程序 clear;x0=input('请输入初值x0:'); ep=input('请输入容许误差:'); flag=1;while flag==1x1=x0-fu(x0)/dfu(x0); if abs(x1-x0)<ep flag=0; end x0=x1; endfprintf('方程的一个近似解为：%f\n',x0); 寻找最大δ值的程序： cleareps=input('请输入搜索精度：'); ep=input('请输入容许误差：'); flag=1; k=0; x0=0; while flag==1 sigma=k*eps; x0=sigma; k=k+1; m=0; flag1=1;while flag1==1 && m<=10^3 x1=x0-fu(x0)/dfu(x0);if abs(x1-x0)<epm=m+1; x0=x1; endif flag1==1||abs(x0)>=ep flag=0; end endfprintf('最大的sigma 值为：%f\n',sigma);2．求下列方程的非零根5130.6651()ln 05130.665114000.0918x x f x x +⎛⎫=-= ⎪-⨯⎝⎭解：Matlab 程序为：（1）主程序 clear clc format long x0=765; N=100;errorlim=10^(-5); x=x0-f(x0)/subs(df(),x0); n=1; while n<Nx=x0-f(x0)/subs(df(),x0); if abs(x-x0)>errorlim n=n+1; else break ; end x0=x; enddisp(['迭代次数: n=',num2str(n)])disp(['所求非零根: 正根x1=',num2str(x),' 负根x2=',num2str(-x)]) （2）子函数 非线性函数f function y=f(x)y=log((513+*x)/*x))-x/(1400*; end（3）子函数 非线性函数的一阶导数df function y=df() syms x1y=log((513+*x1)/*x1))-x1/(1400*; y=diff(y); end运行结果如下：所求非零根: 正根x1= 负根x2=大作业 四试编写MATLAB 函数实现Newton 插值，要求能输出插值多项式. 对函数21()14f x x =+在区间[-5，5]上实现10次多项式插值.分析：（1）输出插值多项式。
北航数值分析大作业三
一、题目：关于x, y, t, u, v, w 的下列方程组0.5cos 2.670.5sin 1.070.5cos 3.740.5sin 0.79t u v w x t u v w y t u v w x t u v w y +++-=⎧⎪+++-=⎪⎨+++-=⎪⎪+++-=⎩1、试用数值方法求出f(x, y)在区域 {(,)|00.8,0.5 1.5}D x y x y =≤≤≤≤上的一个近似表达式,0(,)kr s rsr s p x y cx y ==∑要求(,)p x y 一最小的k 值达到以下的精度10202700((,)(,))10i j i j i j f x y p x y σ-===-≤∑∑其中，0.08,0.50.05i j x i y j ==+。
2、计算****(,),(,)i j i j f x y p x y （i = 1, 2, …，8；j = 1, 2,…，5）的值，以观察(,)p x y 逼近(,)f x y 的效果，其中，*i x =0.1i , *j y =0.5+0.2j 。
说明：1、用迭代方法求解非线性方程组时，要求近似解向量()k x 满足()(1)()12||||/||||10k k k x x x --∞∞-≤2、作二元插值时，要使用分片二次代数插值。
3、要由程序自动确定最小的k 值。
4、打印以下内容：●算法的设计方案。
●全部源程序（要求注明主程序和每个子程序的功能）。
●数表：,,i j x y (,)i j f x y （i = 0,1,2,…,10；j = 0,1,2,…,20）。
●选择过程的,k σ值。
●达到精度要求时的,k σ值以及(,)p x y 中的系数rs c （r = 0,1,…,k;s = 0,1,…,k ）。
●数表：**,,i j x y ****(,),(,)i j i j f x y p x y （i = 1, 2, ...,8；j = 1, 2, (5)。
北航研究生数值分析试题
∗⎞ ⎟的 A1 ⎠
矩阵。
三、（12 分）试用高斯列主元素法求解线性方程组
⎡ 1 3 −2 −4 ⎤ ⎡ x1 ⎤ ⎡3 ⎤ ⎢ 2 6 −7 −10 ⎥ ⎢ x ⎥ ⎢ −2 ⎥ ⎢ ⎥⎢ 2⎥ = ⎢ ⎥ ⎢ −1 −1 5 9 ⎥ ⎢ x3 ⎥ ⎢14 ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥ ⎢ x4 ⎦ ⎥ ⎣ −6 ⎦ ⎣ −3 −5 0 15 ⎦ ⎣ 四、（12 分）利用矩阵 A 的三角分解 A = LU 求解下列方程组 ⎛ 1 2 1 ⎞ ⎛ x1 ⎞ ⎛ 0 ⎞ ⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎜ 2 2 3 ⎟ ⎜ x2 ⎟ = ⎜ 3 ⎟ ⎜ −1 −3 0 ⎟ ⎜ x ⎟ ⎜ 2 ⎟ ⎝ ⎠⎝ 3 ⎠ ⎝ ⎠
第一章
1、近似数 x = 0.231 关于真值 x = 0.229 有（ （1）1；（2）2；（3）3；（4）4。
∗
绪论
一、选择题（四个选项中仅有一项符合题目要求，每小题 3 分，共计 15 分) ）位有效数字。
2、取 3 ≈ 1.732 计算 x = ( 3 − 1) ，下列方法中哪种最好？（
4
）
Ax
∞和
A ∞ 的值分别为（
）
3
(1) 8 , 8 ；
(2) 8 , 7 ；
(3) 8 , 6 ；
(4) 7 , 7 。
5 、若解线性代数方程组的 Gauss 部分选主元方法第二步得到的系数矩阵的第三列向量为
(2
6 3 2 −5 4 2 ) ，则第三步主行是（
T
） (4) 第 6 行。
(1) 第 2 行；
1 − cos x , sin x
x ≠ 0且 x << 1 ；
（2）
1 1− x , − 1+ 2x 1+ x
北航数值分析报告大作业第三题(fortran)
“数值分析“计算实习大作业第三题——SY1415215孔维鹏一、计算说明1、将x i=0.08i，y j=0.5+0.05j分别代入方程组（A.3）得到关于t，u，v，w的的方程组，调用离散牛顿迭代子函数求出与x i，y j对应的t i，u j。
2、调用分片二次代数插值子函数在点(t i,u j)处插值得到z(x i,y j)=f(x i,y j)，得到数表(x i,y j,f(x i,y j))。
3、对于k=1,2,3,4⋯，分别调用最小二乘拟合子函数计算系数矩阵c rs及误差σ，直到满足精度，即求得最小的k值及系数矩阵c rs。
4、将x i∗=0.1i，y j∗=0.5+0.2j分别代入方程组（A.3）得到关于t∗，u∗，v∗，w∗的的方程组，调用离散牛顿迭代子函数求出与x i∗，y j∗对应的t i∗，u j∗，调用分片二次代数插值子函数在点(t i∗,u j∗)处插值得到z∗(x i∗,y j∗)=f(x i∗,y j∗)；调用步骤3中求得的系数矩阵c rs求得p(x i∗,y j∗)，打印数表(x i∗,y j∗,f(x i∗,y j∗),p(x i∗,y j∗))。
二、源程序（FORTRAN）PROGRAM SY1415215DIMENSION X(11),Y(21),T(6),U(6),Z(6,6),UX(11,21),TY(11,21),FXY(11,21),C(6,6) DIMENSION X1(8),Y1(5),FXY1(8,5),PXY1(8,5),UX1(8,5),TY1(8,5)REAL(8) X,Y,T,U,Z,FXY,UX,TY,C,E,X1,Y1,FXY1,PXY1,UX1,TY1OPEN (1,FILE='第三题计算结果.TXT')DO I=1,11X(I)=0.08*(I-1)ENDDODO I=1,21Y(I)=0.5+0.05*(I-1)ENDDO!*****求解非线性方程组，得到z=f（t，u）的函数*******DO I=1,11DO J=1,21CALL DISNEWTON_NONLINEAR(X(I),Y(J),UX(I,J),TY(I,J)) ENDDOENDDO!*************分片二次插值得到z=f（x，y）***********DO I=1,11DO J=1,21CALL INTERPOLATION(UX(I,J),TY(I,J),FXY(I,J))ENDDOENDDOWRITE (1,"('数表(x,y,f(x,y)):')")WRITE (1,"(3X,'X',7X,'Y',10X,'F(X,Y)')")DO I=1,11DO J=1,21WRITE(1,'(1X,F5.2,2X,F5.3,2X,E20.13)') X(I),Y(J),FXY(I,J) ENDDOWRITE (1,"('')")ENDDO!***********最小二乘拟合得到P（x，y）**************N=11M=21WRITE (1,'(" ","K和σ分别为：")')DO K=1,20CALL LSFITTING(X,Y,FXY,C,N,M,K,K,E)WRITE (1,'(I3,2X,E20.13)') K-1,EIF(E<10E-7) EXITENDDOWRITE(1,'(" ")')WRITE(1,'("系数矩阵Crs（按行）为：")')DO I=1,KDO J=1,KWRITE (1,'(E20.13,2X,\)') C(I,J)ENDDOWRITE (1,"('')")WRITE (*,"('')")ENDDODO I=1,8X1(I)=0.1*IENDDODO J=1,5Y1(J)=0.5+0.2*JENDDODO I=1,8DO J=1,5CALL DISNEWTON_NONLINEAR(X1(I),Y1(J),UX1(I,J),TY1(I,J))ENDDOENDDODO I=1,8DO J=1,5CALL INTERPOLATION(UX1(I,J),TY1(I,J),FXY1(I,J))ENDDOENDDOPXY1=0DO I=1,8DO J=1,5DO II=1,KDO JJ=1,KPXY1(I,J)=PXY1(I,J)+C(II,JJ)*(X1(I)**(II-1))*(Y1(J)**(JJ-1)) ENDDOENDDOENDDOENDDOWRITE(1,'(" ")')WRITE(1,'("数表（x,y,f(x,y),p(x,y)):")')WRITE(1,"(2X,'X',6X,'Y',12X,'F(X,Y)',14X,'P(X,Y)')")DO I=1,8DO J=1,5WRITE(1,'(F5.3,2X,F5.3,2X,E20.13,2X,E20.13)') X1(I),Y1(J),FXY1(I,J),PXY1(I,J) ENDDOWRITE (1,"('')")ENDDOCLOSE (1)END!***********用离散牛顿法求解非线性方程组****************SUBROUTINE DISNEWTON_NONLINEAR(X1,Y1,U,T)PARAMETER (N=4)REAL EPS !EPS为迭代精度，M为最大迭代次数DIMENSION X(N),H(N),Y(N),JA(N,N),E(N),XK(N)REAL(8) JA,X,H,Y,E,XK,U,T,V,W,X1,Y1,E1,E2F1(T,U,V,W)=0.5*COS(T)+U+V+W-X1-2.67F2(T,U,V,W)=T+0.5*SIN(U)+V+W-Y1-1.07F3(T,U,V,W)=0.5*T+U+COS(V)+W-X1-3.74F4(T,U,V,W)=T+0.5*U+V+SIN(W)-Y1-0.79EPS=10E-12M=100X=1.0DO K=1,MH=1!计算Y=F（x）Y(1)=F1(X(1),X(2),X(3),X(4))Y(2)=F2(X(1),X(2),X(3),X(4))Y(3)=F3(X(1),X(2),X(3),X(4))Y(4)=F4(X(1),X(2),X(3),X(4))!计算JA（N,N）E=0.0DO I=1,NDO J=1,NDO JJ=1,NIF(JJ==J) THENE(JJ)=X(JJ)+H(JJ)ELSEE(JJ)=X(JJ)ENDIFENDDOIF(I==1) THENJA(I,J)=(F1(E(1),E(2),E(3),E(4))-F1(X(1),X(2),X(3),X(4)))/H(J) ELSEIF(I==2) THENJA(I,J)=(F2(E(1),E(2),E(3),E(4))-F2(X(1),X(2),X(3),X(4)))/H(J) ELSEIF(I==3) THENJA(I,J)=(F3(E(1),E(2),E(3),E(4))-F3(X(1),X(2),X(3),X(4)))/H(J) ELSEIF(I==4) THENJA(I,J)=(F4(E(1),E(2),E(3),E(4))-F4(X(1),X(2),X(3),X(4)))/H(J) ENDIFENDDOENDDO!求解线性方程组CALL GAUSS(JA,XK,-Y,N)!判断精度CALL NORM(XK,N,E1)CALL NORM(X,N,E2)IF(E1/E2<=EPS) THENT=X(1)U=X(2)EXITELSEDO I=1,NX(I)=X(I)+XK(I)ENDDOENDIFENDDORETURNEND!**********列主元高斯消去法求解线性方程组********* SUBROUTINE GAUSS(A,X,B,N)DIMENSION A(N,N),B(N),X(N),T(N,N),TB(N)REAL M(N,N)REAL(8) A,B,X,T!消元过程DO K=1,N-1TA=A(K,K)TL=KDO L=K+1,NIF ((A(L,K)>TA).OR.(A(L,K)==TA)) THENTA=A(L,K)TL=LDO J=K,NT(K,J)=A(K,J)A(K,J)=A(TL,J)A(TL,J)=T(K,J)ENDDOTB(K)=B(K)B(K)=B(TL)B(TL)=TB(K)ENDIFENDDODO I=K+1,NM(I,K)=A(I,K)/A(K,K)A(I,K)=0DO J=K+1,NA(I,J)=A(I,J)-M(I,K)*A(K,J)ENDDOB(I)=B(I)-M(I,K)*B(K)ENDDOENDDO!回代过程X(N)=B(N)/A(N,N)DO K=N-1,1,-1S=0.0DO J=K+1,NS=S+A(K,J)*X(J)ENDDOX(K)=(B(K)-S)/A(K,K)ENDDORETURNEND!***********求向量的无穷数************ SUBROUTINE NORM(X,N,A)DIMENSION X(N)REAL(8) X,AA=ABS(X(1))DO I=2,NIF(ABS(X(I))>ABS(X(I-1))) THENA=ABS(X(I))ENDIFENDDORETURNEND!**************分片二次代数插值************** SUBROUTINE INTERPOLATION(U,V,W) PARAMETER (N=6,M=6)DIMENSION X(N),Y(M),Z(M,N),LK(3),LR(3)REAL(8) X,Y,Z,H,TREAL(8) U,V,W,LK,LR !U,V分别为插值点处的坐标，W为插值结果INTEGER R!**********************数据赋值********************** DATA Y/0.0,0.2,0.4,0.6,0.8,1.0/DATA X/0.0,0.4,0.8,1.2,1.6,2.0/DATA Z/-0.5,-0.42,-0.18,0.22,0.78,1.5,&&-0.34,-0.5,-0.5,-0.34,-0.02,0.46,&&0.14,-0.26,-0.5,-0.58,-0.5,-0.26,&&0.94,0.3,-0.18,-0.5,-0.66,-0.66,&&2.06,1.18,0.46,-0.1,-0.5,-0.74,&&3.5,2.38,1.42,0.62,-0.02,-0.5/H=0.4T=0.2!******************计算K,R*************************IF(U<=X(2)+H/2) THENK=2ELSEIF(U>X(N-1)-H/2) THENK=N-1ELSEDO I=3,N-2IF((U>X(I)-H/2).AND.(U<=X(I)+H/2)) THENK=IENDIFENDDOENDIFIF(V<=Y(2)+T/2) THENR=2ELSEIF(V>Y(M-1)-T/2) THENR=M-1ELSEDO J=3,M-2IF((V>Y(J)-T/2).AND.(V<=Y(J)+T/2)) THENR=JENDIFENDDOENDIFI=KJ=RLK(1)=(U-X(I))*(U-X(I+1))/(X(I-1)-X(I))/(X(I-1)-X(I+1))LK(2)=(U-X(I-1))*(U-X(I+1))/(X(I)-X(I-1))/(X(I)-X(I+1)) LK(3)=(U-X(I))*(U-X(I-1))/(X(I+1)-X(I))/(X(I+1)-X(I-1)) LR(1)=(V-Y(J))*(V-Y(J+1))/(Y(J-1)-Y(J))/(Y(J-1)-Y(J+1)) LR(2)=(V-Y(J-1))*(V-Y(J+1))/(Y(J)-Y(J-1))/(Y(J)-Y(J+1)) LR(3)=(V-Y(J))*(V-Y(J-1))/(Y(J+1)-Y(J))/(Y(J+1)-Y(J-1))W=0DO K=1,3DO R=1,3W=W+LK(K)*LR(R)*Z(J+R-2,I+K-2)ENDDOENDDORETURNEND!*******************最小二乘拟合子函数************** SUBROUTINE LSFITTING(X,Y,Z,A,N,M,P,Q,DT1)INTEGER P,QDIMENSION X(N),Y(M),Z(N,M),A(P,Q)DIMENSION APX(20),APY(20),BX(20),BY(20),U(20,20),V(20,M) DIMENSION T(20),T1(20),T2(20)REAL(8) X,Y,Z,A,DT1DO I=1,PDO J=1,QA(I,J)=0.0ENDDOENDDOIF(P>N) P=NIF(P>20) P=20IF(Q>M) Q=MIF(Q>20) Q=20XX=0YY=0D1=NAPX(1)=0.0DO I=1,NAPX(1)=APX(1)+X(I)ENDDOAPX(1)=APX(1)/D1DO J=1,MV(1,J)=0.0DO I=1,NV(1,J)=V(1,J)+Z(I,J)ENDDOV(1,J)=V(1,J)/D1ENDDOIF(P>1) THEND2=0.0APX(2)=0.0DO I=1,NG=X(I)-APX(1)D2=D2+G*GAPX(2)=APX(2)+(X(I)-XX)*G*G ENDDOAPX(2)=APX(2)/D2BX(2)=D2/D1DO J=1,MV(2,J)=0.0DO I=1,NG=X(I)-APX(1)V(2,J)=V(2,J)+Z(I,J)*GENDDOV(2,J)=V(2,J)/D2ENDDOD1=D2ENDIFDO K=3,PD2=0.0APX(K)=0.0DO J=1,MV(K,J)=0.0ENDDODO I=1,NG1=1.0G2=X(I)-APX(1)DO J=3,KG=(X(I)-APX(J-1))*G2-BX(J-1)*G1 G1=G2G2=GENDDOD2=D2+G*GAPX(K)=APX(K)+X(I)*G*GDO J=1,MV(K,J)=V(K,J)+Z(I,J)*GENDDOENDDODO J=1,MV(K,J)=V(K,J)/D2ENDDOAPX(K)=APX(K)/D2BX(K)=D2/D1D1=D2ENDDOD1=MAPY(1)=0.0DO I=1,MAPY(1)=APY(1)+Y(I)ENDDOAPY(1)=APY(1)/D1DO J=1,PU(J,1)=0.0DO I=1,MU(J,1)=U(J,1)+V(J,I) ENDDOU(J,1)=U(J,1)/D1ENDDOIF(Q>1)THEND2=0.0APY(2)=0.0DO I=1,MG=Y(I)-APY(1)D2=D2+G*GAPY(2)=APY(2)+(Y(I))*G*G ENDDOAPY(2)=APY(2)/D2BY(2)=D2/D1DO J=1,PU(J,2)=0.0DO I=1,MG=Y(I)-APY(1)U(J,2)=U(J,2)+V(J,I)*GENDDOU(J,2)=U(J,2)/D2ENDDOD1=D2ENDIFDO K=3,QD2=0.0APY(K)=0.0DO J=1,PU(J,K)=0.0ENDDODO I=1,MG1=1.0G2=Y(I)-APY(1)DO J=3,KG=(Y(I)-APY(J-1))*G2-BY(J-1)*G1 G1=G2G2=GENDDOD2=D2+G*GAPY(K)=APY(K)+Y(I)*G*GDO J=1,PU(J,K)=U(J,K)+V(J,I)*GENDDOENDDODO J=1,PU(J,K)=U(J,K)/D2ENDDOAPY(K)=APY(K)/D2BY(K)=D2/D1D1=D2ENDDOV(1,1)=1.0V(2,1)=-APY(1)V(2,2)=1.0DO I=1,PDO J=1,QA(I,J)=0.0ENDDOENDDODO I=3,QV(I,I)=V(I-1,I-1)V(I,I-1)=-APY(I-1)*V(I-1,I-1)+V(I-1,I-2)IF(I>=4) THENDO K=I-2,2,-1V(I,K)=-APY(I-1)*V(I-1,K)+V(I-1,K-1)-BY(I-1)*V(I-2,K) ENDDOENDIFV(I,1)=-APY(I-1)*V(I-1,1)-BY(I-1)*V(I-2,1)ENDDODO I=1,PIF(I==1) THENT(1)=1.0T1(1)=1.0ELSEIF(I==2) THENT(1)=-APX(1)T(2)=1.0T2(1)=T(1)T2(2)=T(2)ELSET(I)=T2(I-1)T(I-1)=-APX(I-1)*T2(I-1)+T2(I-2)IF(I>=4) THENDO K=I-2,2,-1T(K)=-APX(I-1)*T2(K)+T2(K-1)-BX(I-1)*T1(K)ENDDOENDIFT(1)=-APX(I-1)*T2(1)-BX(I-1)*T1(1)T2(I)=T(I)DO K=I-1,1,-1T1(K)=T2(K)T2(K)=T(K)ENDDOENDIFDO J=1,QDO K=I,1,-1DO L=J,1,-1A(K,L)=A(K,L)+U(I,J)*T(K)*V(J,L)ENDDOENDDOENDDOENDDODT1=0.0DO I=1,NX1=X(I)DO J=1,MY1=Y(J)X2=1.0DD=0.0DO K=1,PG=A(K,Q)DO KK=Q-1,1,-1G=G*Y1+A(K,KK)ENDDOG=G*X2DD=DD+GX2=X2*X1ENDDODT=DD-Z(I,J)DT1=DT1+DT*DTENDDOENDDORETURNEND三、计算结果数表(x,y,f(x,y)):X Y UX TY F(X,Y) 0.00 0.500 1.345 0.243 0.17E+000.00 0.550 1.322 0.269 0.66E+000.00 0.600 1.299 0.295 0.35E+000.00 0.650 1.277 0.322 0.94E+000.00 0.700 1.255 0.350 0.30E-020.00 0.750 1.235 0.377 -0.87E-010.00 0.800 1.215 0.406 -0.58E+000.00 0.850 1.196 0.434 -0.72E+000.00 0.900 1.177 0.463 -0.54E+000.00 0.950 1.159 0.492 -0.86E+000.00 1.050 1.125 0.550 -0.74E+00 0.00 1.100 1.109 0.580 -0.06E+00 0.00 1.150 1.093 0.609 -0.00E+00 0.00 1.200 1.079 0.639 -0.18E+00 0.00 1.250 1.064 0.669 -0.52E+00 0.00 1.300 1.050 0.699 -0.19E+00 0.00 1.350 1.037 0.729 -0.48E+00 0.00 1.400 1.024 0.759 -0.68E+00 0.00 1.450 1.011 0.790 -0.52E+00 0.00 1.500 1.000 0.820 -0.29E+000.08 0.500 1.415 0.228 0.67E+00 0.08 0.550 1.391 0.253 0.08E+00 0.08 0.600 1.368 0.279 0.02E+00 0.08 0.650 1.346 0.306 0.47E+00 0.08 0.700 1.325 0.333 0.57E+00 0.08 0.750 1.304 0.360 0.48E-01 0.08 0.800 1.284 0.388 -0.73E-01 0.08 0.850 1.265 0.416 -0.16E+00 0.08 0.900 1.246 0.444 -0.29E+00 0.08 0.950 1.229 0.473 -0.36E+00 0.08 1.000 1.211 0.502 -0.08E+00 0.08 1.050 1.194 0.531 -0.29E+00 0.08 1.100 1.178 0.560 -0.78E+00 0.08 1.150 1.163 0.589 -0.93E+00 0.08 1.200 1.148 0.619 -0.44E+00 0.08 1.250 1.133 0.649 -0.92E+00 0.08 1.300 1.119 0.679 -0.71E+000.08 1.400 1.093 0.739 -0.37E+00 0.08 1.450 1.080 0.769 -0.83E+00 0.08 1.500 1.068 0.799 -0.92E+000.16 0.500 1.483 0.214 0.31E+00 0.16 0.550 1.460 0.239 0.64E+00 0.16 0.600 1.437 0.264 0.91E+00 0.16 0.650 1.414 0.290 0.06E+00 0.16 0.700 1.393 0.316 0.70E+00 0.16 0.750 1.372 0.343 0.59E+00 0.16 0.800 1.352 0.370 0.12E+00 0.16 0.850 1.333 0.398 0.77E-02 0.16 0.900 1.315 0.426 -0.83E-01 0.16 0.950 1.297 0.454 -0.58E+00 0.16 1.000 1.279 0.483 -0.20E+00 0.16 1.050 1.262 0.512 -0.11E+00 0.16 1.100 1.246 0.541 -0.74E+00 0.16 1.150 1.231 0.570 -0.09E+00 0.16 1.200 1.216 0.600 -0.59E+00 0.16 1.250 1.201 0.629 -0.66E+00 0.16 1.300 1.187 0.659 -0.71E+00 0.16 1.350 1.174 0.689 -0.32E+00 0.16 1.400 1.161 0.718 -0.56E+00 0.16 1.450 1.148 0.748 -0.31E+00 0.16 1.500 1.136 0.778 -0.75E+000.24 0.500 1.551 0.201 0.66E+01 0.24 0.550 1.527 0.225 0.03E+000.24 0.650 1.482 0.275 0.64E+00 0.24 0.700 1.460 0.301 0.47E+00 0.24 0.750 1.439 0.327 0.34E+00 0.24 0.800 1.419 0.354 0.24E+00 0.24 0.850 1.400 0.381 0.69E+00 0.24 0.900 1.381 0.409 0.04E-01 0.24 0.950 1.363 0.437 -0.42E-01 0.24 1.000 1.346 0.465 -0.06E+00 0.24 1.050 1.329 0.494 -0.59E+00 0.24 1.100 1.313 0.523 -0.83E+00 0.24 1.150 1.297 0.552 -0.15E+00 0.24 1.200 1.282 0.581 -0.19E+00 0.24 1.250 1.267 0.610 -0.84E+00 0.24 1.300 1.253 0.640 -0.66E+00 0.24 1.350 1.240 0.669 -0.30E+00 0.24 1.400 1.227 0.699 -0.86E+00 0.24 1.450 1.214 0.729 -0.84E+00 0.24 1.500 1.202 0.759 -0.77E+000.32 0.500 1.617 0.188 0.28E+01 0.32 0.550 1.593 0.212 0.49E+01 0.32 0.600 1.570 0.236 0.68E+00 0.32 0.650 1.547 0.261 0.75E+00 0.32 0.700 1.526 0.286 0.60E+00 0.32 0.750 1.505 0.312 0.77E+00 0.32 0.800 1.485 0.339 0.05E+00 0.32 0.850 1.466 0.365 0.99E+00 0.32 0.900 1.447 0.393 0.27E+000.32 1.000 1.411 0.448 -0.01E-02 0.32 1.050 1.395 0.477 -0.41E-01 0.32 1.100 1.378 0.505 -0.18E+00 0.32 1.150 1.363 0.534 -0.25E+00 0.32 1.200 1.347 0.563 -0.29E+00 0.32 1.250 1.333 0.592 -0.90E+00 0.32 1.300 1.319 0.621 -0.00E+00 0.32 1.350 1.305 0.650 -0.40E+00 0.32 1.400 1.292 0.680 -0.54E+00 0.32 1.450 1.279 0.710 -0.79E+00 0.32 1.500 1.267 0.739 -0.91E+000.40 0.500 1.681 0.177 0.91E+01 0.40 0.550 1.658 0.199 0.00E+01 0.40 0.600 1.634 0.223 0.83E+01 0.40 0.650 1.612 0.247 0.02E+01 0.40 0.700 1.591 0.272 0.94E+00 0.40 0.750 1.570 0.298 0.49E+00 0.40 0.800 1.550 0.324 0.94E+00 0.40 0.850 1.530 0.350 0.40E+00 0.40 0.900 1.512 0.377 0.33E+00 0.40 0.950 1.493 0.405 0.99E+00 0.40 1.000 1.476 0.432 0.68E+00 0.40 1.050 1.459 0.460 0.08E-01 0.40 1.100 1.443 0.488 -0.84E-01 0.40 1.150 1.427 0.517 -0.98E+00 0.40 1.200 1.412 0.545 -0.27E+00 0.40 1.250 1.397 0.574 -0.06E+000.40 1.350 1.369 0.632 -0.66E+00 0.40 1.400 1.356 0.662 -0.37E+00 0.40 1.450 1.343 0.691 -0.43E+00 0.40 1.500 1.331 0.721 -0.12E+000.48 0.500 1.745 0.166 0.69E+01 0.48 0.550 1.721 0.188 0.02E+01 0.48 0.600 1.698 0.211 0.74E+01 0.48 0.650 1.676 0.235 0.40E+01 0.48 0.700 1.654 0.259 0.23E+01 0.48 0.750 1.634 0.284 0.56E+00 0.48 0.800 1.613 0.310 0.28E+00 0.48 0.850 1.594 0.336 0.49E+00 0.48 0.900 1.575 0.363 0.31E+00 0.48 0.950 1.557 0.390 0.66E+00 0.48 1.000 1.539 0.417 0.30E+00 0.48 1.050 1.522 0.444 0.34E+00 0.48 1.100 1.506 0.472 0.07E-01 0.48 1.150 1.490 0.500 -0.62E-01 0.48 1.200 1.475 0.529 -0.45E+00 0.48 1.250 1.460 0.557 -0.86E+00 0.48 1.300 1.446 0.586 -0.39E+00 0.48 1.350 1.432 0.615 -0.22E+00 0.48 1.400 1.419 0.644 -0.67E+00 0.48 1.450 1.406 0.674 -0.55E+00 0.48 1.500 1.394 0.703 -0.14E+000.56 0.500 1.808 0.156 0.48E+010.56 0.600 1.761 0.200 0.10E+01 0.56 0.650 1.739 0.223 0.68E+01 0.56 0.700 1.717 0.247 0.94E+01 0.56 0.750 1.696 0.272 0.33E+01 0.56 0.800 1.676 0.297 0.11E+00 0.56 0.850 1.657 0.323 0.63E+00 0.56 0.900 1.638 0.349 0.97E+00 0.56 0.950 1.620 0.375 0.52E+00 0.56 1.000 1.602 0.402 0.56E+00 0.56 1.050 1.585 0.429 0.47E+00 0.56 1.100 1.568 0.457 0.20E+00 0.56 1.150 1.552 0.485 0.13E+00 0.56 1.200 1.537 0.513 0.09E-01 0.56 1.250 1.522 0.541 -0.47E-01 0.56 1.300 1.508 0.570 -0.99E+00 0.56 1.350 1.494 0.599 -0.82E+00 0.56 1.400 1.481 0.627 -0.26E+00 0.56 1.450 1.468 0.657 -0.71E+00 0.56 1.500 1.455 0.686 -0.98E+000.64 0.500 1.870 0.147 0.74E+01 0.64 0.550 1.846 0.168 0.10E+01 0.64 0.600 1.823 0.190 0.54E+01 0.64 0.650 1.801 0.213 0.42E+01 0.64 0.700 1.779 0.236 0.56E+01 0.64 0.750 1.758 0.260 0.03E+01 0.64 0.800 1.738 0.285 0.42E+01 0.64 0.850 1.718 0.310 0.41E+010.64 0.950 1.681 0.362 0.36E+00 0.64 1.000 1.664 0.388 0.18E+00 0.64 1.050 1.646 0.415 0.28E+00 0.64 1.100 1.630 0.443 0.07E+00 0.64 1.150 1.614 0.470 0.66E+00 0.64 1.200 1.598 0.498 0.09E+00 0.64 1.250 1.584 0.526 0.50E-01 0.64 1.300 1.569 0.554 -0.88E-01 0.64 1.350 1.555 0.583 -0.76E+00 0.64 1.400 1.542 0.611 -0.66E+00 0.64 1.450 1.529 0.640 -0.33E+00 0.64 1.500 1.516 0.669 -0.56E+000.72 0.500 1.931 0.139 0.94E+01 0.72 0.550 1.907 0.159 0.84E+01 0.72 0.600 1.884 0.181 0.36E+01 0.72 0.650 1.862 0.203 0.40E+01 0.72 0.700 1.840 0.226 0.47E+01 0.72 0.750 1.819 0.249 0.56E+01 0.72 0.800 1.799 0.273 0.19E+01 0.72 0.850 1.779 0.298 0.37E+01 0.72 0.900 1.760 0.323 0.86E+01 0.72 0.950 1.742 0.349 0.76E+00 0.72 1.000 1.724 0.375 0.24E+00 0.72 1.050 1.707 0.402 0.55E+00 0.72 1.100 1.691 0.429 0.97E+00 0.72 1.150 1.675 0.456 0.27E+00 0.72 1.200 1.659 0.484 0.31E+000.72 1.250 1.644 0.511 0.51E+00 0.72 1.300 1.630 0.539 0.49E+00 0.72 1.350 1.616 0.568 0.72E-02 0.72 1.400 1.602 0.596 -0.69E-01 0.72 1.450 1.589 0.625 -0.67E+00 0.72 1.500 1.576 0.653 -0.20E+000.80 0.500 1.992 0.131 0.31E+01 0.80 0.550 1.968 0.151 0.44E+01 0.80 0.600 1.945 0.172 0.41E+01 0.80 0.650 1.922 0.193 0.45E+01 0.80 0.700 1.900 0.216 0.00E+01 0.80 0.750 1.879 0.239 0.10E+01 0.80 0.800 1.859 0.263 0.16E+01 0.80 0.850 1.840 0.287 0.52E+01 0.80 0.900 1.821 0.312 0.02E+01 0.80 0.950 1.802 0.337 0.38E+01 0.80 1.000 1.784 0.363 0.89E+01 0.80 1.050 1.767 0.389 0.28E+00 0.80 1.100 1.751 0.416 0.09E+00 0.80 1.150 1.734 0.443 0.23E+00 0.80 1.200 1.719 0.470 0.93E+00 0.80 1.250 1.704 0.498 0.15E+00 0.80 1.300 1.689 0.525 0.86E+00 0.80 1.350 1.675 0.553 0.64E+00 0.80 1.400 1.662 0.582 0.74E-01 0.80 1.450 1.649 0.610 -0.37E-01 0.80 1.500 1.636 0.638 -0.81E+00K和σ分别为：0 0.93E+031 0.61E+012 0.92E-023 0.53E-034 0.16E-055 0.77E-07系数矩阵Crs（按行）为：0.00E+01 -0.83E+01 0.56E+00 0.97E+00 -0.03E+00 0.70E-010.91E+01 -0.99E+00 -0.96E+01 0.17E+01 -0.66E+00 0.10E-010.77E+00 0.42E+01 -0.10E+00 -0.81E+00 0.81E+00 -0.62E-01-0.25E+00 -0.21E+00 0.97E+00 -0.18E+00 0.49E+00 -0.63E-010.34E+00 -0.56E+00 0.69E-01 0.51E+00 -0.77E-01 0.27E-01-0.94E-01 0.94E+00 -0.58E+00 0.69E-01 -0.50E-01 0.53E-02数表（x,y,f(x,y),p(x,y)):X Y F(X,Y) P(X,Y)0.100 0.700 0.58E+00 0.05E+000.100 1.100 -0.66E+00 -0.26E+00 0.100 1.300 -0.68E+00 -0.31E+00 0.100 1.500 -0.52E+00 -0.49E+000.200 0.700 0.54E+00 0.19E+00 0.200 0.900 -0.63E-01 -0.65E-01 0.200 1.100 -0.90E+00 -0.90E+00 0.200 1.300 -0.84E+00 -0.90E+00 0.200 1.500 -0.03E+00 -0.04E+000.300 0.700 0.82E+00 0.09E+00 0.300 0.900 0.48E+00 0.11E+00 0.300 1.100 -0.63E+00 -0.88E+00 0.300 1.300 -0.72E+00 -0.96E+00 0.300 1.500 -0.34E+00 -0.84E+000.400 0.700 0.79E+00 0.89E+00 0.400 0.900 0.56E+00 0.63E+00 0.400 1.100 -0.83E-01 -0.04E-01 0.400 1.300 -0.72E+00 -0.71E+00 0.400 1.500 -0.85E+00 -0.07E+000.500 0.700 0.56E+01 0.92E+01 0.500 0.900 0.51E+00 0.23E+00 0.500 1.100 0.59E+00 0.27E+00 0.500 1.300 -0.53E+00 -0.11E+00 0.500 1.500 -0.67E+00 -0.33E+000.600 0.900 0.14E+00 0.75E+00 0.600 1.100 0.19E+00 0.32E+00 0.600 1.300 -0.70E-01 -0.82E-01 0.600 1.500 -0.08E+00 -0.75E+000.700 0.700 0.89E+01 0.29E+01 0.700 0.900 0.91E+01 0.11E+01 0.700 1.100 0.60E+00 0.97E+00 0.700 1.300 0.22E-01 0.06E-01 0.700 1.500 -0.53E+00 -0.80E+000.800 0.700 0.09E+01 0.06E+01 0.800 0.900 0.32E+01 0.50E+01 0.800 1.100 0.03E+00 0.79E+00 0.800 1.300 0.25E+00 0.50E+00 0.800 1.500 -0.14E+00 -0.28E+00。
北航数值分析全部三次大作业
北航数值分析全部三次大作业第一次大作业是关于解线性方程组的数值方法。
我们被要求实现各种常用的线性方程组求解算法，例如高斯消元法、LU分解法和迭代法等。
我首先学习了这些算法的原理和实现方法，并借助Python编程语言编写了这些算法的代码。
在实验中，我们使用了不同规模和条件的线性方程组进行测试，并比较了不同算法的性能和精度。
通过这个作业，我深入了解了线性方程组求解的原理和方法，提高了我的编程和数值计算能力。
第二次大作业是关于数值积分的方法。
数值积分是数值分析中的重要内容，它可以用于计算曲线的长度、函数的面积以及求解微分方程等问题。
在这个作业中，我们需要实现不同的数值积分算法，例如矩形法、梯形法和辛普森法等。
我学习了这些算法的原理和实现方法，并使用Python编写了它们的代码。
在实验中，我们计算了不同函数的积分值，并对比了不同算法的精度和效率。
通过这个作业，我深入了解了数值积分的原理和方法，提高了我的编程和数学建模能力。
第三次大作业是关于常微分方程的数值解法。
常微分方程是数值分析中的核心内容之一，它可以用于描述众多物理、化学和生物现象。
在这个作业中，我们需要实现不同的常微分方程求解算法，例如欧拉法、龙格-库塔法和Adams法等。
我学习了这些算法的原理和实现方法，并使用Python编写了它们的代码。
在实验中，我们解决了一些具体的常微分方程问题，并比较了不同算法的精度和效率。
通过这个作业，我深入了解了常微分方程的原理和方法，提高了我的编程和问题求解能力。
总的来说，北航数值分析课程的三次大作业非常有挑战性，但也非常有意义。
通过这些作业，我在数值计算和编程方面得到了很大的提升，也更加深入地了解了数值分析的理论和方法。
虽然这些作业需要大量的时间和精力，但我相信这些努力将会对我未来的学习和工作产生积极的影响。
BUAA数值分析大作业三
北京航空航天大学2020届研究生《数值分析》实验作业第九题院系：xx学院学号：姓名：2020年11月Q9:方程组A.4一、 算法设计方案（一）总体思路1.题目要求∑∑===k i kj s r rsy x cy x p 00),(对f(x, y) 进行拟合，可选用乘积型最小二乘拟合。
),(i i y x 与),(i i y x f 的数表由方程组与表A-1得到。
2.),(**j i y x f 与1使用相同方法求得，),(**j i y x p 由计算得出的p(x,y)直接带入),(**j i y x 求得。
1. ),(i i y x 与),(i i y x f 的数表的获得对区域D ={ （x,y)|1≤x ≤1.24,1.0≤y ≤1.16}上的f (x , y )值可通过xi=1+0.008i ，yj=1+0.008j ，得到),(i i y x 共31×21组。
将每组带入A4方程组，即可获得五个二元函数组，通过简单牛顿迭代法求解这五个二元数组可获得z1~z5有关x,y 的表达式。
再将),(i i y x 分别带入z1~z5表达式即可获得f(x,y)值。
2.乘积型最小二乘曲面拟合2.1使用乘积型最小二乘拟合，根据k 值不用，有基函数矩阵如下：⎪⎪⎪⎭⎫ ⎝⎛=k i i k x x x x B 0000 ， ⎪⎪⎪⎭⎫ ⎝⎛=k j jk y y y y G 0000数表矩阵如下：⎪⎪⎪⎭⎫⎝⎛=),(),(),(),(0000j i i j y x f y x f y x f y x f U记C=[rs c ]，则系数rs c 的表达式矩阵为：11-)(-=G G UG B B B C T TT ）（通过求解如下线性方程，即可得到系数矩阵C 。
UG B G G C B B T T T =)()(2.2计算),(),,(****j i j i y x p y x f (i =1,2,…,31 ; j =1,2,…,21) 的值),(**j i y x f 的计算与),(j i y x f 相同。
北航数值分析大作业3(学硕)
《数值分析》作业三院系：机械学院学号：SY1307145姓名：龙安林2013年11 月24 日1. 算法设计1) 开始；2) 计算数组[][]0.08,0.050.5,0,1,2,,10;0,1,2,,20x i i y j j i j ==+=⋯=⋯（）； 3) 将点[][],0,1,2,,10;0,1,2,,20x i y j i j =⋯=⋯（）,（）带入非线性方程组： 0.5cos 2.670.5sin 1.070.5cos 3.740.5sin 0.79t u v w x t u v w y t u v w x t u v w y +++-=⎧⎪+++-=⎪⎨+++-=⎪⎪+++-=⎩ 得出相应的点,t u （）； 4) 选择拉格朗日插值法，将,t u （）作为中间变量，在题目所给出的二维数表中进行二次代数插值，得到[][],)(z f x i y j =；5) 输出数表：[][][][]()()0,1,2,,10;0,1,2,,20,,,x i y j f x i y j i j =⋯=⋯； 6) 令k=0；7) 以()()(),,,0,1,r r r s x x y y r s ϕψ===…,k 为拟合基函数，将上述数表作为拟合条件，对于给定的k 值，得到矩阵B 、G 、U ；8) 令-1-1(),()T T T A B B B U C AG G G ==，用选主元的LU 分解法分别计算矩阵A 和C 的各列，最后得到系数矩阵C ；9) 以公式：()()()00,k ki j rs r i s j s r p x y C x y ϕψ===∑∑计算每个点的拟合值；10) 利用公式：()()()2102000,,i j i j i j f x y p x y σ===-∑∑计算拟合误差，当σ≤10-7时，循环结束，否则k=k+1，转（6）；11) 令[][]()**0.10.50.2 1,2,81,2,5x i i y j j i j ==+=⋯=⋯；，；，；12) 计算()()()******,,,,,i j i j i jf x y p x y delta x y ，输出数表，观察逼近效果； 13) 结束。
北航研究生数值分析大作业三
数值分析—计算实习作业三学院：17系专业：精密仪器及机械姓名：张大军学号：DY1417114一、程序设计方案程序设计方案流程图如图1所示。
（注：由本人独立完成，并且有几处算法很巧妙，但同时也有许多不足，可以优化和模块化，由于时间原因只实现了调试通过）图1.程序设计方案流程图二、程序源代码#include <iostream.h>#include <iomanip.h>#include <math.h>#include<stdio.h>#include <conio.h>#define M 10000#define N 4#define E 1.0e-12int zuixiaci;static double c[9][9];static double bijin[8][5];int main(){double X[N]={0,0,0,1};double T[11][21],U[11][21],xianshi[11][21];double diertX[N];double F[N];double f[N][N];double Max1=0,Max2=0;int k,i,j,t,tt=0,yao=0;void qiuF(double * X,double *F,int i,int j);void qiuF2(double *X,double *F,int i,int j);void qiuf(double * X,double (*f)[N]);void qiudiertX(double (*a)[N],double*b,double*X); double gouzaohs(double t,double u); void solve_C(double (*U)[21]); void pp(double (*U)[21],int k);for(i=0;i<11;i++)for(j=0;j<21;j++){for(k=0;k<M;k++){qiuF(X,F,i,j);qiuf(X,f);qiudiertX(f,F,diertX);for(t=0;t<N;t++){X[t]=X[t]+diertX[t];}Max1=0,Max2=0;for(t=0;t<N;t++){if(Max1<fabs(X[t]))Max1=fabs(X[t]);if(Max2<fabs(diertX[t]))Max2=fabs(diertX[t]);}if((Max2/Max1)<=E){k=M;yao=1;T[i][j]=X[0];U[i][j]=X[1];xianshi[i][j]=gouzaohs(X[0],X[1]);cout<<setiosflags(ios::scientific)<<setprecision(12);cout<<setprecision(2)<<"("<<setw(5)<<0.08*i<<","<<setw(5)<< (0.5+0.05*j)<<",";cout<<setprecision(12)<<setw(21)<<xianshi[i][j]<<") ";if(tt==3){tt=0;cout<<'\n';cout<<'\n';}else{tt++;}}}if(yao==0)cout<<"迭代不成功"<<endl; yao=0;}cout<<endl;solve_C(xianshi);pp(xianshi,zuixiaci);tt=0;for(i=1;i<9;i++)for(j=1;j<6;j++){for(k=0;k<M;k++){qiuF2(X,F,i,j);qiuf(X,f);qiudiertX(f,F,diertX);for(t=0;t<N;t++){X[t]=X[t]+diertX[t];}Max1=0,Max2=0;for(t=0;t<N;t++){if(Max1<fabs(X[t]))Max1=fabs(X[t]);if(Max2<fabs(diertX[t]))Max2=fabs(diertX[t]);}if((Max2/Max1)<=E){k=M;yao=1;xianshi[i-1][j-1]=gouzaohs(X[0],X[1]);cout<<setiosflags(ios::scientific)<<setprecision(12);cout<<setprecision(2)<<"("<<setw(5)<<0.1*i<<","<<setw(5)<<( 0.5+0.2*j)<<",";cout<<setprecision(12)<<setw(21)<<xianshi[i-1][j-1]<<","<<set w(21)<<bijin[i-1][j-1]<<") ";if(tt==2){tt=0;cout<<'\n';}else{tt++;}}}if(yao==0)cout<<"迭代不成功"<<endl;yao=0;}cout<<endl;return 1;}void qiuF(double *X,double *F,int i,int j){F[0]=-(0.5*cos(X[0])+X[1]+X[2]+X[3]-0.08*i-2.67);F[1]=-(X[0]+0.5*sin(X[1])+X[2]+X[3]-(0.5+0.05*j)-1.07);F[2]=-(0.5*X[0]+X[1]+cos(X[2])+X[3]-0.08*i-3.74);F[3]=-(X[0]+0.5*X[1]+X[2]+sin(X[3])-(0.5+0.05*j)-0.79); }void qiuF2(double *X,double *F,int i,int j){F[0]=-(0.5*cos(X[0])+X[1]+X[2]+X[3]-0.1*i-2.67);F[1]=-(X[0]+0.5*sin(X[1])+X[2]+X[3]-(0.5+0.2*j)-1.07);F[2]=-(0.5*X[0]+X[1]+cos(X[2])+X[3]-0.1*i-3.74);F[3]=-(X[0]+0.5*X[1]+X[2]+sin(X[3])-(0.5+0.2*j)-0.79); }void qiuf(double *X,double (*f)[N]){f[0][0]=-0.5*sin(X[0]);f[0][1]=1;f[0][2]=1;f[0][3]=1;f[1][0]=1;f[1][1]=0.5*cos(X[1]);f[1][2]=1;f[1][3]=1;f[2][0]=0.5;f[2][1]=1;f[2][2]=-sin(X[2]);f[2][3]=1;f[3][0]=1;f[3][1]=0.5;f[3][2]=1;f[3][3]=cos(X[3]);}//求解关于变化X的线性方程组void qiudiertX(double (*a)[N],double*b,double*X) {double H[N][N]={0},l[N]={0};double B;double sum;int i,j,m,k,z;for(k=0;k<N-1;k++){for(j=k;j<N;j++){l[j]=a[k][j];}z=k;for(m=k;m<N;m++){if(fabs(a[z][k])<fabs(a[m][k]))z=m;}for(j=k;j<N;j++){a[k][j]=a[z][j];a[z][j]=l[j];}B=b[k];b[k]=b[z];b[z]=B;for(i=k+1;i<N;i++){H[i][k]=a[i][k]/a[k][k];for(j=k+1;j<N;j++)a[i][j]=a[i][j]-H[i][k]*a[k][j];b[i]=b[i]-H[i][k]*b[k];}}if(a[N-1][N-1]==0){cout<<"算法失效，停止计算"<<endl; }else{X[N-1]=b[N-1]/a[N-1][N-1];for(k=N-2;k>=0;k--){sum=0;for(j=k+1;j<N;j++){sum=sum+a[k][j]*X[j];}X[k]=(b[k]-sum)/a[k][k];}}}//作二元差值，使用分片二次代数插值double gouzaohs(double t,double u){double T[6]={0,0.2,0.4,0.6,0.8,1},U[6]={0,0.4,0.8,1.2,1.6,2};double Z[6][6]={-0.5,-0.34,0.14,0.94,2.06,3.5,-0.42,-0.5,-0.26,0.3,1.18,2.38,-0.18,-0.5,-0.5,-0.18,0.46,1.42,0.22,-0.34,-0.58,-0.5,-0.1,0.62,0.78,-0.02,-0.5,-0.66,-0.5,-0.02,1.5,0.46,-0.26,-0.66,-0.74,-0.5};double g=0,sum=0,sum1=1,sum2=1;int i=0,j=0,k=0,r=0,kk=0,rr=0;for(i=1;(i<6)&&(T[i]-0.1<t);i++){}for(j=1;(j<6)&&(U[j]-0.2<u);j++){}if(i==1)i=2;if(i==6)i=5;if(j==1)j=2;if(j==6)j=5;sum=0;for(k=i-2;k<i+1;k++)for(r=j-2;r<j+1;r++){sum1=1;sum2=1;for(kk=i-2;kk<i+1;kk++){if(k!=kk){sum1=sum1*(t-T[kk])/(T[k]-T[kk]);}}for(rr=j-2;rr<j+1;rr++){if(r!=rr){sum2=sum2*(u-U[rr])/(U[r]-U[rr]);}}sum=sum+sum1*sum2*Z[k][r];}g=sum;return g;}//求r*s阶矩阵A与s*t阶矩阵B的乘积矩阵Cvoid Multi(double *a, double *b, double *c, int la, int lb, int lc, int r, int s, int t){int i, j, k;for (i=0; i<r; i++)for (j=0; j<t; j++){*(c+i*lc+j)=0;for (k=0; k<s; k++)*(c+i*lc+j)+=*(a+i*la+k)*(*(b+k*lb+j));}}//求n阶方阵A的逆矩阵Bdouble Inverse(double *a, double *b, int la, int lb, int n){int i, j, k;double temp;for(i=0; i<n; i++)for(j=0; j<n; j++)if (i==j)*(b+i*lb+j)=1;else*(b+i*lb+j)=0;for (k=0; k<n; k++){j=k;for (i=k+1; i<n; i++)if (fabs(*(a+i*la+k))>fabs(*(a+j*la+k))) j=i;if (j!=k)for (i=0; i<n; i++){temp=*(a+j*la+i);*(a+j*la+i)=*(a+k*la+i);*(a+k*la+i)=temp;temp=*(b+j*lb+i);*(b+j*lb+i)=*(b+k*lb+i);*(b+k*lb+i)=temp;}if (*(a+k*la+k)==0)return 0;if ((temp=*(a+k*la+k))!=1)for (i=0; i<n; i++){*(a+k*la+i)/=temp;*(b+k*lb+i)/=temp;}for (i=0; i<n; i++)if ((*(a+i*la+k)!=0) && (i!=k)){temp=*(a+i*la+k);for (j=0; j<n; j++){*(a+i*la+j)-=temp*(*(a+k*la+j));*(b+i*lb+j)-=temp*(*(b+k*lb+j));}}}return 0;}void solve_C(double (*U)[21]){int i,j,r,s,k;double t1[21][21], t2[21][21], t3[21][21],d[9][9],e[9][9];double B[11][9], B_T[9][11], G[21][9], G_T[9][21],P[11][21];double temp, FangCha;for(i=0;i<9;i++){for(j=0;j<11;j++){B[j][i]=pow(0.08*j,i);B_T[i][j]=pow(0.08*j,i);}for(j=0;j<21;j++){G[j][i]=pow(0.5+0.05*j,i);G_T[i][j]=pow(0.5+0.05*j,i);}}for (k=0; k<9; k++){FangCha=0;Multi(B_T[0], B[0], t1[0], 11, 9, 21, k+1, 11, k+1);Inverse(t1[0], c[0], 21, 9, k+1);Multi(e[0], c[0], d[0], 9, 9, 9, k+1, k+1, k+1);Multi(c[0], B_T[0], t1[0], 9, 11, 21, k+1, k+1, 11);Multi(t1[0], U[0], t2[0], 21, 21, 21, k+1, 11, 21);Multi(G_T[0], G[0], t1[0], 21, 9, 21, k+1, 21, k+1);Inverse(t1[0], c[0], 21, 9, k+1);Multi(G[0], c[0], t3[0], 9, 9, 21, 21, k+1, k+1);Multi(t2[0], t3[0], c[0], 21, 21, 9, k+1, 21, k+1);for(i=0;i<11;i++)for(j=0;j<21;j++){temp=0;for(r=0;r<k+1;r++)for(s=0;s<k+1;s++)temp+=c[r][s]*B[i][r]*G[j][s];P[i][j]=temp;FangCha+=(U[i][j]-temp)*(U[i][j]-temp);}cout<<"k="<<setw(5)<<k<<";"<<setw(5)<<"Sigma="<<FangCha<<" ;\n"<<'\n';if(FangCha<=1.0e-7){zuixiaci=k;cout<<"达到精度要求时: k="<<setw(5)<<k<<";"<<setw(5)<<"Sigma="<<FangCha<<";\n";cout<<" 系数c(r,s)如下：\n";for(i=0;i<k+1;i++){for(j=0;j<k+1;j++){cout<<"C("<<i<<","<<j<<")="<<setw(21)<<c[i][j]<<"; ";}cout<<endl<<'\n';}cout<<endl;return;}}cout<<"经过8次拟合没有达到所需精度；"<<endl;//最高可拟合10次return;}void pp(double (*U)[21],int k){int i,j,r,s;double B[8][9],G[5][9],temp;for(i=0;i<k+1;i++){for(j=0;j<8;j++){B[j][i]=pow(0.1*(j+1),i);}for(j=0;j<5;j++){G[j][i]=pow(0.5+0.2*(j+1),i);}}for(i=0;i<8;i++)for(j=0;j<5;j++){temp=0;for(r=0;r<k+1;r++)for(s=0;s<k+1;s++)temp+=c[r][s]*B[i][r]*G[j][s];bijin[i][j]=temp;}}三、程序运行结果显示程序运行结果显示如图2。
北航《数值分析》习题
北航《数值分析》习题习题一1. 下列各数都是经过四舍五入得到的近似数，试指出它们有几位有效数字以及它们的绝对误差限、相对误差限。
（1）；（2）；（3）；（4）；（5）；（6）；（7）；1. （1）5，，；（2）2，，；（3）4，，；（4）5，，；（5）1，，；（6）2，，；（7）6，，2. 为使下列各数的近似值的相对误差限不超过，问各近似值分别应取几位有效数字？2. ；；3. 设均为第1题所给数据，估计下列各近似数的误差限。
（1）；（2）；（3）3. （1）；（2）；（3）4. 计算，取，利用下列等价表达式计算，哪一个的结果最好？为什么？（1）；（2）；（3）（4）4. 第（3）个结果最好。
5. 序列满足递推关系式若（三位有效数字），计算时误差有多大？这个计算过程稳定吗？5. 不稳定。
从计算到时，误差约为6. 求方程的两个根，使其至少具有四位有效数字（要求利用。
6. ，7. 某生产部门生产的一件产品需用七个零件，而这七个零件的质量取决于零件参数的标定值，它们的参数允许有一定的误差：若每一零件的标定值取做区间中点，在生产过程中每一零件的参数都有可能产生误差。
由此将零件分成不同的等级：A，B，C三等，等级由标定值的相对误差限表示，A等为1%，B等为5%，C等为10%。
试确定三个等级的零件分别满足的区间。
8. 将一个八位二进制数(10111101)2转换成十进制数时，可以用公式：（1）用多项式求值的秦九韶方法求C的值；（2）写出将任意一个八位二进制数(b1b2b3b4b5b6b7b8)2转化为十进制数的算法。
9. 利用等式变换使下列表达式的计算结果比较精确。
（1）；（2）（3）；（4）9. （1）；（2）；（3）；（4）10. 设，求证：（1）（2）利用（1）中的公式正向递推计算时误差增大；反向递推时误差函数减小。
习题二1. 判断下列方程有几个实根，并求出其隔根区间。
（1）；（2）（3）；（4）1. （1），，；（2）；（3），，；（4）为根。
北航数分大作业三
一、算法的设计方案1、对于已给出的非线性方程组，其解集可采用牛顿迭代法进行求解。
在每次迭代过程中，将x ，y 的值固定，如此便可得到一组关于t ，u ，v ，w 的解。
因此可以建立一组（x ，y ）和（t ，u ）一一对应的关系。
2、采用分片二次插值对题目中所给出的z ，t ，u 二维数表进行处理。
于是在 0≤t ≤1, 0≤u ≤2 的矩形区域就建立了 z 与(t,u)的一一对应关系。
其中选择（m ，n ）满足,2322m i m h h t t t m -<≤+≤≤，,2322n j n u u u n ττ-<≤+≤≤。
3、对i x i *08.0= 10,...,2,1,0=i ，j y j 05.05.0+= 20,...,2,1,0=j 。
分别使用前两步算法，可得到一组2111)(],[⨯=j i y x f j i z 的数表。
4、采用最小二乘拟合，设∑∑===k r k s s r rs y x cy x p 00)(，m=10 n=20，M=N=K 。
插值基函数10,...,1,0,)(==i x x r i i r φ k r ,...,1,0=，ks j y y sj j s ,...,1,020,...,1,0,)(===ψ。
、)1()1()1()1(][][+⨯++⨯+==k n sj k m r i y G x B U 即为上面所求的Z[11][21]。
为避免计算过程中出现矩阵求逆，将U B B B A T T 1)(-=改为U B A B B T T =)(，再利用高斯消去法以)(B B T作为系数矩阵，U B T 的每一列作为非线性部分，分别解出A 的每一列。
在将1)(-=G G AG A T 改为AG G G C T =)(，然后利用高斯消元法以)(G G T 作为系数矩阵，AG 的每一行作为非线性部分，分解出C 的每一行。
如此便得到了最小二乘拟合的系数矩阵C 。
北航数值分析B自测题1-3章
D.2
10.若f (x)在区间[a, b]内单调有限，二分k 次后区间记为[ak , bk ],且每次 取xk+1 = ak 近似代替精确解x∗ ,则最小的绝对误差限为（ ） A.|xk+1 − x∗ | ≤ C.|xk+1 − x∗ | ≤ 二、 填空题 1.设π 的近似数π ∗ 有4位有效数字，则其相对误差限为 √ 2. x∗ 的相对误差约是x∗ 的相对误差的 倍。
(2)a=6,b=2; (3)a=2,b=3; (4)a=-1,b=2.
3.设矩阵A ∈ Rn∗n , Q ∈ Rn∗n ,且QT Q=E，则下列关系式不成立的是（ ） (1) A
2
= AQ ; (2) QA
2 F
= A
F
; (3) Qx
2
= x , 其中x ∈ Rn ;
2
(4)cond∞ ( ) = cond∞ ( Q).
8
3 −1 4 1 4.设矩阵A=−1 2 −2 ，x=−1 ，则 Ax 2 −3 −2 1 （ ） (1)8,8; (2)8,7; (3)8,6; (4)7,7.
A. 7 3
B.− 7 3
7 C. 6
D.− 7 6
8.追赶法适用于求解（ ）线性方程组. A.上三角 B.下三角 C.对角 D.三对角 二.填空 题：
2 1 0 1.设A=1 2 a，为使A 可分解为A=LLT ，其中L是 对角元素为正的下 0 a 2 三角矩阵 范围是 ，则a的取值 .
1 (s∗ ) s∗
≈
27 110×80
= 0.31% = max 1 ≤ j ≤ n{1, 2, 4} = 4;
北航数值分析大作业题目三
《数值分析》第三次大作业一、算法的设计方案： （一）、总体方案设计：（1）解非线性方程组。
将给定的(,)i i x y 当作已知量代入题目给定的非线性方程组，求得与(,)i i x y 相对应的数组t[i][j],u[i][j]。
（2）分片二次代数插值。
通过分片二次代数插值运算，得到与数组t[11][21],u[11][21]]对应的数组z[11][21]，得到二元函数z=(,)i i f x y 。
（3）曲面拟合。
利用x[i],y[j],z[11][21]建立二维函数表，再根据精度的要求选择适当k 值，并得到曲面拟合的系数矩阵C[r][s]。
（4）观察和(,)i i p x y 的逼近效果。
观察逼近效果只需要重复上面（1）和（2）的过程，得到与新的插值节点(,)i i x y 对应的(,)i i f x y ，再与对应的(,)i i p x y 比较即可，这里求解(,)i i p x y 可以直接使用（3）中的C[r][s]和k 。
（二）具体算法设计：（1）解非线性方程组牛顿法解方程组()0F x =的解*x ，可采用如下算法： 1）在*x 附近选取(0)xD ∈，给定精度水平0ε>和最大迭代次数M 。
2）对于0,1,k M =L 执行① 计算()()k F x和()()k F x '。
② 求解关于()k x∆的线性方程组()()()()()k k k F x x F x '∆=-③ 若()()k k x x ε∞∞∆≤，则取*()k x x ≈，并停止计算；否则转④。
④ 计算(1)()()k k k xx x +=+∆。
⑤ 若k M <，则继续，否则，输出M 次迭代不成功的信息，并停止计算。
（2）分片双二次插值给定已知数表以及需要插值的节点，进行分片二次插值的算法：设已知数表中的点为： 00(0,1,,)(0,1,,)i j x x ih i n y y j j m τ=+=⎧⎪⎨=+=⎪⎩L L ，需要插值的节点为(,)x y 。
北航数值分析大作业题目三
1、 算法的设计方案： （一）、总体方案设计：
（1）解非线性方程组。将给定的当作已知量代入题目给定的非线性方程组，求得与相对应 的数组t[i][j],u[i][j]。 （2）分片二次代数插值。通过分片二次代数插值运算，得到与数组t[11][21],u[11][21]] 对应的数组z[11][21]，得到二元函数z=。 （3）曲面拟合。利用x[i],y[j],z[11][21]建立二维函数表，再根据精度的要求选择适当k 值，并得到曲面拟合的系数矩阵C[r][s]。 （4）观察和的逼近效果。观察逼近效果只需要重复上面（1）和（2）的过程，得到与新的 插值节点对应的，再与对应的比较即可，这里求解可以直接使用（3）中的C[r][s]和k。
{ temp=0; for(l=k+1;l<=3;l++) {temp=temp+dF[k][l]*dx[l]/dF[k][k];} dx[k]=-F[k]/dF[k][k]-temp; } temp=0; for(l=0;l<=3;l++) /*求解矩阵范数，用无穷范数*/ { if(temp<fabs(dx[l])) temp=fabs(dx[l]); } fx=temp; temp=0; for(l=0;l<=3;l++) { if(temp<fabs(X[l])) temp=fabs(X[l]); } fX=temp; if(fabs(fx/fX)<Epsilon1) /*判断是否成立*/ { t[i][j]=X[0]; u[i][j]=X[1]; goto loop4;} else { for(l=0;l<=3;l++) {X[l]=X[l]+dx[l];} n=n+1; goto loop3;} } loop3:{if(n<M) /*判断是否超出规定迭代次数*/ goto loop1; else printf("迭代不成功\n"); goto loop4; } loop4:{continue;} } } } void fpeccz(double t[11][21],double u[11][21])/*分片二次代数插值子程序*/ { int s[11][21],r[11][21]; int i,j,i1,j1,m; double z0[6][6]={{-0.5,-0.34,0.14,0.94,2.06,3.5}, {-0.42,-0.5,-0.26,0.3,1.18,2.38}, {-0.18,-0.5,-0.5,-0.18,0.46,1.42},
北航数值分析A大作业3
一、算法设计方案1、解非线性方程组将各拟合节点（x i ，y j ）分别带入非线性方程组，求出与(,)i i x y 相对应的数组te[i][j]，ue[i][j]，求解非线性方程组选择Newton 迭代法，迭代过程中需要求解线性方程组，选择选主元的Doolittle 分解法。
2、二元二次分偏插值对数表z(t,u)进行分片二次代数插值，求得对应（t ij ，u ij ）处的值，即为),(j i y x f 的值。
根据给定的数表，可将整个插值区域分成 16 个小 的区域，故先判断(t i j , u ij ) 所在，的区域，再作此区域的插值，计算 z ij ，相应的Lagrange 形式的插值多项式为：°112211(,)()()(,)m n krkrk m r n p t u l t l u f t u ++=-=-=∑∑其中11()m wk w m k ww kt t l t t t +=-≠-=-∏ (k=m-1, m, m+1) °11()n wr w n r ww ry y l u y y +=-≠-=-∏ (r=n-1, n, n+1)3、曲面拟合从k=1开始逐渐增大k 的值，使用最小二乘法曲面拟合法对z=f(x,y)进行拟合，当710-<σ时结束计算。
拟合基函数φr (x)ψs (y)选择为φr (x)=x r ，ψs (y)=y s 。
拟合系数矩阵c 通过连续两次解线性方程组求得。
[]rsc *=C ，11()()T T T --=C B B B UG G G其中011101011[()]1kk r i k x x x x x x x ϕ⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎣⎦B L LM M M M L ，0011101011[()]1k k s j k y y y y G y y y ψ⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎣⎦L LM M M M L [(,)]i j f x y =U4、观察比较计算)5,,2,1,8,,2,1)(,(),,(****⋅⋅⋅=⋅⋅⋅=j i y x p y x f j i j i 的值并输出结果，以观察),(y x p 逼近),(y x f 的效果。
北航数值分析复习试题
数值分析一、单项选择题（共20分,每小题2分）1-110=11=12=，则Lagranage 二次插值多项式为（ ） A.2(121)(144)(100)(144)(100)(121)()101112(100121)(100144)(121100)(121144)(144121)(144100)x x x x x x L x ------=++------ B ．2(121)(144)(100)(144)(100)(121)()111012(100121)(100144)(121100)(121144)(144121)(144100)x x x x x x L x ------=++------ C ．2(121)(144)(100)(144)(100)(121)()121110(100121)(100144)(121100)(121144)(144121)(144100)x x x x x x L x ------=++------ D ．2(121)(144)(100)(144)(100)(121)()101211(100121)(100144)(121100)(121144)(144121)(144100)x x x x x x L x ------=++------ 1-210=11=12=，用Lagranage为（ ）精确到小数点后4位。
A.9.7227 B ．11.7227 C ．10.7227 D ．13.72271-3、已知(1 2 3 4)TX =，则向量X 的21, , Xx x ∞的值分别是：（ ）B. -9，212,7C. 4,5,6D. 9,4,71-4、设 2121A --⎛⎫= ⎪⎝⎭，则21,, , F A A A x ∞的值分别为（ ）A. 4B. -9，C.4,5,6D. 9,4,71-5、设节点00 (=0,1,2,...,n), (0),k x x kh k x x th t =+=+>则Newton 向前插值公式为（ ）A. 100010()()!k k nn k j f N x th f t j k -==∆+=+-∑∏ B. 110()()!k k nn n n n k j f N x th f t j k -==∆+=+-∑∏ C. 100010()()!k k nn k j f N x th f t j k -==∇+=+-∑∏ D. 110()()!k k nn n n n k j f N x th f t j k -==∇+=+-∑∏1-6、方程组⎪⎪⎩⎪⎪⎨⎧=+++=+++=+++=+++47401815622189622315694962424321432143214321x x x x x x x x x x x x x x x x 进行直接三角分解法得到的L 矩阵为（ ）A. 1211213321B.1613216241C.16332202102 D.12147165511-7、对方程组的系数矩阵123412312341346262414535x x x x x x x x x x x x x x ++-=⎧⎪++=-⎨++-=⎪--+=-⎩进行Crout 分解法得到的U 矩阵为（ ）A.1111363111261131-- B. 1111363111569171--C.111136611151091371-- D. 111166311223611121--1-8、1、已知642()1f x x x x =+-+，2, 2 (0,1,2,...)k x kh h k =+==，则[2,6,10,14,18,22,26,30]f =（ ）A ．5！B ．4！C ．0D ．11-9、1、已知64()f x x x =+，2, 2 (0,1,2,...)k x kh h k =+==，则[2,4,6,8,10,12,14]f =（ ）A ．5！B ．4！C ．0D ．11-10、复合Cotes 求积公式, 复合梯形求积公式和复合Simpson 求积公式的收敛阶分别为（ ） A ．5，1，3 B ．4，2 ，6 C ．6，2，4 D ．以上都不对1-11、对线性方程组1231231232211221x x x x x x x x x +-=⎧⎪++=⎨++=⎪⎩，若用Jocabi 迭代法和G-S 迭代法求解，则（ ）A.Jocabi 迭代法收敛和G-S 迭代法发散B. Jocabi 迭代法和G-S 迭代法均发散C. Jocabi 迭代法和G-S 迭代法均收敛D. Jocabi 迭代法发散和G-S 迭代法收敛1-12、对线性方程组1231213918 293x x x x x x x --=⎧⎪-+=⎨-+=⎪⎩，若用Jocabi 迭代法和G-S 迭代法求解（ ），则 B.Jocabi 迭代法收敛和G-S 迭代法发散 A. Jocabi 迭代法和G-S 迭代法均发散C. Jocabi 迭代法和G-S 迭代法均收敛D. Jocabi 迭代法发散和G-S 迭代法收敛1-13、设线性方程组为1231213918 293x x x x x x x --=⎧⎪-+=⎨-+=⎪⎩，则Jocabi 迭代格式和G-S 迭代格式分别为（ ），则(Ⅰ) 2311(1)()()1(1)()2(1)()311799917881899k k k k k k k x x x x x x x +++⎧=++⎪⎪⎪=+⎨⎪⎪=+⎪⎩(Ⅱ) 2311(1)()()1(1)(1)2(1)(1)311799917881899k k k k k k k x x x x x x x +++++⎧=++⎪⎪⎪=+⎨⎪⎪=+⎪⎩A.(Ⅰ)和(Ⅱ)B. (Ⅱ)和(Ⅰ)C.(Ⅰ)和(Ⅰ)D. (Ⅱ)和(Ⅱ)1-14、已知*x 是()f x 的 (2)m m ≥重根，则求重根的修正Newton 公式为（ ）1(). ()k k k k f x A x x mf x +=-' 10(). ()k k k f x B x x mf x +=-'111(). ()()()k k k k k k k f x C x x x x f x f x +--=--- 111()(). ()()k k k k k k k f x f x D x x f x x x -+--=--1-15、若记(),()k k k k y f x z f y ==，则对迭代格式1()k k x f x -=使用Aitken 加速后得到的新迭代迭代格式为（ ）21(()). (())2()k k k k k k k f x x A x x f f x f x x +-=--+21(()). ()(())2()k k k k k k k f x x B x f x f f x f x x +-=--+21(). 2k k k k k k k z y C x z z y x +-=--+ 21((())()). (())(())2()k k k k k k k f f x f x D x f f x f f x f x x +-=--+1-16、将积分区间[a,b]n 等分，分点为kh a x k +=，k=0,1,2,3,4....n,其中nab h -=，则复合梯形公式为( )A. ])()(4)([211∑-=++n k k b f x f a f hB.])()(2)([211∑-=++n k k b f x f a f hC.)]()(4)(2)([6102111b f x f x f a f hn k k n k k +++∑∑-=+-=D.)]()(4)(2)([6112110b f x f x f a f hn k k n k k +++∑∑-=+-=二、填空题（共20分,每空2分）2-1、根据数值方法的稳定性与算法设计原则在连加运算中要防止 ，在减法运算中要避免 ，在除法运算中要避免，在乘法运算中要避免 。
北航研究生数值分析期末模拟历年考试
北航研究⽣数值分析期末模拟历年考试数值分析模拟试卷1⼀、填空（共30分，每空3分） 1 设-=1511A ，则A 地谱半径=)(a ρ______，A 地条件数)(1A cond =________. 2设,2,1,0,,53)(2==+=k kh x x x f k ，则],,[21++n n n x x x f ＝________,],,[321+++n n n n x x x x f ，＝________.3 设≤≤-++≤≤+=21,1210,)(2323x cx bx x x x x x S ，是以0，1，2为节点地三次样条函数，则b=________,c=________.4设∞=0)]([k k x q 是区间［0，1］上权函数为x x =)(ρ地最⾼项系数为1地正交多项式族，其中1)(0=x q ，则=1)(dx x xq k________，=)(2x q ________.5设=11001a a a a A ，当∈a ________时，必有分解式，其中L 为下三⾓阵，当其对⾓线元素)3,2,1(=i L ii 满⾜条件________时，这种分解是唯⼀地.⼆、（14分）设49,1,41,)(21023====x x x x x f , （1）试求)(x f 在]49,41[上地三次Hermite 插值多项式)(x H 使满⾜2,1,0),()(==i x f x H i i ，)()(11x f x H '='.（2）写出余项)()()(x H x f x R -=地表达式.三、（14分）设有解⽅程0cos 2312=+-x x 地迭代公式为n n x x cos 3241+=+，（1）证明R x ∈?0均有?∞→=x x n x lim （?x 为⽅程地根）；（2）取40=x ，⽤此迭代法求⽅程根地近似值，误差不超过，列出各次迭代值；（3）此迭代地收敛阶是多少？证明你地结论.四、(16分) 试确定常数A ，B ，C 和，使得数值积分公式有尽可能⾼地代数精度. 试问所得地数值积分公式代数精度是多少？它是否为Gauss 型地？五、（15分）设有常微分⽅程地初值问题??=='00)(),(y x y y x f y ，试⽤Taylor 展开原理构造形如)()(11011--++++=n n n n n f f h y y y ββα地⽅法，使其具有⼆阶精度，并推导其局部截断误差主项.六、（15分）已知⽅程组b Ax ＝，其中==21,13.021b A ，（1）试讨论⽤Jacobi 迭代法和Gauss-Seidel 迭代法求解此⽅程组地收敛性.（2）若有迭代公式)()()()1(b Ax a x x k k k ++=+，试确定⼀个地取值范围，在这个范围内任取⼀个值均能使该迭代公式收敛.七、（8分）⽅程组，其中，A 是对称地且⾮奇异.设A 有误差，则原⽅程组变化为，其中为解地误差向量，试证明.其中1λ和2λ分别为A 地按模最⼤和最⼩地特征值.数值分析模拟试卷2填空题（每空2分，共30分）1. 近似数231.0=*x 关于真值229.0=x 有____________位有效数字；设)(x f 可微，求⽅程)(x f x =根地⽜顿迭代格式是_______________________________________________；3. 对1)(3++=x x x f ，差商=]3,2,1,0[f _________________；=]4,3,2,1,0[f ________；4. 已知-='-=1223,)3,2(A x ，则=∞||||Ax ________________，=)(1A Cond ______________________ ；⽤⼆分法求⽅程01)(3=-+=x x x f 在区间[0,1]内地根，进⾏⼀步后根所在区间为_________，进⾏⼆步后根所在区间为_________________；求解线性⽅程组=+=+04511532121x x x x 地⾼斯—赛德尔迭代格式为_______________________________________；该迭代格式迭代矩阵地谱半径=)(G ρ_______________；为使两点数值求积公式：-+≈111100)()()(x f x f dx x f ωω具有最⾼地代数精确度，其求积节点应为=0x _____ ,=1x _____,==10ωω__________.8. 求积公式)]2()1([23)(3f f dx x f +≈?是否是插值型地__________，其代数精度为___________.⼆、（12分）(1)设LU A =，其中L 为下三⾓阵，U 为单位上三⾓阵.已知------=2100121001210012A ，求L，U . (2)设A 为66?矩阵，将A 进⾏三⾓分解：LU A =，L 为单位下三⾓阵，U 为上三⾓阵，试写出L 中地元素65l 和U 中地元素56u 地计算公式.三、（12分）设函数)(x f 在区间[0,3]上具有四阶连续导数,试确定⼀个次数不超过3地多项式)(x H ，满⾜3)1()1(,1)2()2(,1)1()1(,0)0()0(='='======f H f H f H f H ，并写出插值余项. （12分）线性⽅程组=+=-22112122b x x b x x ρρ（1）请写出解此⽅程组地赛德尔迭代法地迭代格式，并讨论收敛性. （2）设2=ρ，给定松弛因⼦21=ω，请写出解此⽅程组地SOR ⽅法地迭代格式，并讨论收敛性.五、（7分）改写⽅程042=-+x x为2ln /)4ln(x x -=地形式，问能否⽤迭代法求所给⽅程在[1,2]内地实根？六、（7分）证明解⽅程0)(23=-a x 求3a 地⽜顿迭代法仅为线性收敛. 七、（12分）已知.43,21,41210===x x x （1）推导以这3个点作为求积节点在[0,1]上地插值型求积公式；（2）指明求积公式具有地代数精度；（3）⽤所求公式计算12dx x.⼋、（8分）若i n x x x x x x x x f ),())(()(10---= 互异，求],,,[10p x x x f 地值，这⾥.1+≤n p数值分析模拟试卷3⼀、填空题（每空3分，共30分）1．设1234)(248+++=x x x x f ，则差商=]2,,2,2[810 f ； 2．在⽤松弛法(SOR)解线性⽅程组b Ax =时，若松弛因⼦ω满⾜1|1|≥-ω，则迭代法；3．设,0)(,0)(**≠'=x f x f 要使求*x 地Newton 迭代法⾄少三阶收敛，)(x f 需要满⾜；4. 设)133)(2()(23-+-+=x x x x x f ,⽤Newton 迭代法求21-=x 具有⼆阶收敛地迭代格式为________________ ；求12=x 具有⼆阶收敛地迭代格式为___________________；5．已知?--=1327A ，则=)(A ρ__________,=∞)(A Cond ______ 6. 若1>>x ，改变计算式1lg lg 2--x x =___________________，使计算结果更为精确； 7．过节点())3,2,1,0(,3=i x x i i 地插值多项式为_____________ ； 8. 利⽤抛物(Simpson)公式求212dx x =.⼆、（14分）已知⽅阵=123111122A ，(1) 证明： A 不能被分解成⼀个单位下三⾓阵L 和⼀个上三⾓阵U 地乘积；(2) 给出A 地选主元地Doolittle 分解，并求出排列阵；(3) ⽤上述分解求解⽅程组b Ax =，其中Tb )4,2,5.3(=.三、（12分）设函数)(x f 在区间[0,3]上具有四阶连续导数,试确定⼀个次数不超过3地多项式)(x H ，满⾜40)1()1(,10)1()1(,1)1()1(,0)0()0(=''=''='='-====f H f H f H f H ，并写出插值余项.四、（10分）证明对任意地初值0x ，迭代格式n n x x cos 1=+均收敛于⽅程x x cos =地根，且具有线性收敛速度.五、（12分）在区间[-1,1]上给定函数14)(3+=x x f ，求其在},,1{2x x Span =φ中关于权函数1)(=x ρ地最佳平⽅逼近多项式.（可⽤数据：2123)(,)(,1)(2210-===x x p x x p x p ）六、(12分)(1)试导出切⽐雪夫(Chebyshev)正交多项式])1,1[,,2,1,0)(arccos cos()(-∈==x n x n x T n 地三项递推关系式：=-===-+),2,1()()(2)(,)(,1)(1110 n x T x xT x T x x T x T n n n （2）⽤⾼斯—切⽐雪夫求积公式计算积分dx x x x I ? --=22)2(1，问当节点数n 取何值时，能得到积分地精确值？并计算它.七、（10分）验证对?-+-+=++==++=?+))1(,)1((),(),()(2,13121311hK t y h t x f K thK y th x f K y x f K K K h y y t n n n n n n n n 为2阶格式.参考答案1 ⼀、1．6)(=a ρ，)(1A cond =6.2．],,[21++n n n x x x f =3，],,[321+++n n n n x x x x f ，=0. 3．b =－2，c=3.4．??≠=0,00,21k k ；10356)(22+-=x x x q .5．)3,2,1(0);21,21(=>-∈i l a ii⼆、(1) 25145023345026322514)(23-++-=x x x x H (2) ).49,41(),49()1)(41(169!41)(225∈---=-ξξx x x x R三、（1）32=L ;（2）347.3≈?x ;（3）线性收敛. 四、512,916,910-====αB C A ;求积公式具有5次代数精度，是Gauss 型地. 五、41472110＝－，＝，＝ββα；截断误差主项为)(833n x y h '''. 六、（1）,16.0)(,6.0)(<==G S J B B ρρ因此两种迭代法均收敛.（2）当06.011>>+a 时，该迭代公式收敛.参考答案2 ⼀、1．22．),1,0()()(1 ='-=+n x f x f x x n n n n3．1, 0 4．7,725 5．)43,21(),1,21( 6. 121,2013531)1(1)1(2)(2)1(1??-=-=+++k k k k x x x x 7. 32,3210=-=x x ; 18. 是, 1⼆、(1)---=---=100431000321000211,4510003410002310002U L (2))(;)(4654356532652165155565545643563256215616565u l u l u l u l a u u u l u l u l u l a l +++-=+++-=三、)2()1(!4)()(),2)(1(2)(2)4(--=---=x x x f x R x x x x x H ξ四、(1) ??-=+=+++)1(12)1(2)(21)1(12k k k k x b x x b x ρρ, 1<ρ时收敛(2) ??-+=++=+++)1(1)(22)1(2)(2)(11)1(1214212k k k k k k x x b x x x b x , 收敛五、收敛七、（1）)43(32)21(31)41(32f f f +- （2）2 (3)31 ⼋、110时为时为+=≤n ，p n p参考答案3 ⼀、1．42．发散3．0)(*=''x f4．),1,0()()(1 ='-=+n x f x f x x n n n n ,),1,0()()(31 ='-=+n x f x f x x n n n n5．2608+, 49 6.1lg2-x x7. 3x 8.37 ⼆、(2) 先交换2、3两⾏，交换1、2两⾏，===010001100,5.0003333.06667.00123,15.03333.0016667.0001P U L(3) )5.4,1,5.1('-三、3)4(2)1(!4)()(,)1(9)1(11)(-=-+-+-=x x f x R x x x x x x H ξ五、10512p p +六、1=n ，2π版权申明本⽂部分内容，包括⽂字、图⽚、以及设计等在⽹上搜集整理.版权为个⼈所有This article includes some parts, including text, pictures, and design. Copyright is personal ownership.SixE2yXPq5⽤户可将本⽂地内容或服务⽤于个⼈学习、研究或欣赏，以及其他⾮商业性或⾮盈利性⽤途，但同时应遵守著作权法及其他相关法律地规定，不得侵犯本⽹站及相关权利⼈地合法权利.除此以外，将本⽂任何内容或服务⽤于其他⽤途时，须征得本⼈及相关权利⼈地书⾯许可，并⽀付报酬.6ewMyirQFLUsers may use the contents or services of this article for personal study, research or appreciation, and othernon-commercial or non-profit purposes, but at the same time, they shall abide by the provisions of copyright law and other relevant laws, and shall not infringe upon the legitimate rights of this website and its relevant obligees. In addition, when any content or service of this article is used for other purposes, written permission and remuneration shall be obtained from the person concerned and the relevant obligee.kavU42VRUs转载或引⽤本⽂内容必须是以新闻性或资料性公共免费信息为使⽤⽬地地合理、善意引⽤，不得对本⽂内容原意进⾏曲解、修改，并⾃负版权等法律责任.y6v3ALoS89Reproduction or quotation of the content of this article must be reasonable and good-faith citation for the use of news or informative public free information. It shall not misinterpret or modify the original intention of the content of this article, and shall bear legal liability such as copyright.M2ub6vSTnP。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
{0.22,-0.34,-0.58,-0.5,-0.1,0.62}, {0.78,-0.02,-0.5,-0.66,-0.5,-0.02}, {1.5,0.46,-0.26,-0.66,-0.74,-0.5}}; double t0[6]={0,0.2,0.4,0.6,0.8,1.0}; double u0[6]={0,0.4,0.8,1.2,1.6,2.0}; double temp1,temp2; for(i=0;i<=10;i++) /*选取插值节点*/ { for(j=0;j<=20;j++) { if(t[i][j]<=0.3) s[i][j]=1; else if(t[i][j]>0.7) s[i][j]=4; else { for(m=2;m<=3;m++) { if((t[i][j]>0.2*m-0.1)&&(t[i][j]<=0.2*m+0.1)) { s[i][j]=m; } } } } } for(i=0;i<=10;i++) { for(j=0;j<=20;j++) { f(u[i][j]>1.4) r[i][j]=4; else { for(m=2;m<=3;m++) { if((u[i][j]>0.4*m-0.2)&&(u[i][j]<=0.4*m+0.2)) { r[i][j]=m; } } } } } for(i=0;i<=10;i++) /*插值运算*/ { for(j=0;j<=20;j++) { z[i][j]=0; for(i1=s[i][j]-1;i1<=s[i][j]+1;i1++)
dF[0][0]=-0.5*sin(X[0]); /*求解F'（x）*/ dF[0][1]=1; dF[0][2]=1; dF[0][3]=1; dF[1][0]=1; dF[1][1]=0.5*cos(X[1]); dF[1][2]=1; dF[1][3]=1; dF[2][0]=0.5; dF[2][1]=1; dF[2][2]=-sin(X[2]); dF[2][3]=1; dF[3][0]=1; dF[3][1]=0.5; dF[3][2]=1; dF[3][3]=cos(X[3]); /*高斯选主元消去法求解Δx*/ for(k=0;k<3;k++) { ik=k; for(l=k;l<=3;l++) {if(dF[ik][k]<dF[l][k]) ik=l; } /*选主元*/ temp=0; temp=F[ik]; F[ik]=F[k]; F[k]=temp; for(l=k;l<=3;l++) { temp=0; temp=dF[ik][l]; dF[ik][l]=dF[k][l]; dF[k][l]=temp; } for(l=k+1;l<=3;l++) { m=dF[l][k]/dF[k][k]; F[l]=F[l]-m*F[k]; for(p=k+1;p<=3;p++) {dF[l][p]=dF[l][p]-m*dF[k][p];} } /*消去过程*/ } dx[3]=-F[3]/dF[3][3]; for(k=2;k>=0;k--)
{ for(j1=r[i][j]-1;j1<=r[i][j]+1;j1++) { temp1=1.0; for(m=s[i][j]-1;m<=s[i][j]+1;m++) { if(m!=i1) temp1*=(t[i][j]-t0[m])/(t0[i1]-t0[m]); } temp2=1.0; for(m=r[i][j]-1;m<=r[i][j]+1;m++) { if(m!=j1) temp2*=(u[i][j]-u0[m])/(u0[j1]-u0[m]); } z[i][j]+=temp1*temp2*z0[i1][j1]; } } } } } void qmnh() /*曲面拟合子程序*/ { int i,j,k,m,l,i1,j1,ik; double A[N][21],D[N][21],B[11][N],Bt[N][11],BtB[N][N],BtU[N][21],BtB1[N][N]; double G[21][N],Gt[N][21],GtG[N][N],GtG1[N][N]; double sigma,p[11][21],temp,q; printf("选择过程的k和sigma值:\n"); k=0; sigma=1; /*选取初值，使循环开始*/ /*求解系数矩阵Crs*/ while(sigma>1e-7) { for(i=0;i<=10;i++) { for(j=0;j<=k;j++) { B[i][j]=pow(x[i],j); Bt[j][i]=B[i][j]; } } for(i=0;i<=k;i++) { for(j=0;j<=k;j++) { temp=0; for(l=0;l<=10;l++) { temp+=Bt[i][l]*B[l][j]; } BtB[i][j]=temp; } } for(i=0;i<=k;i++) { for(j=0;j<=20;j++)
《数值分析》第三次大作业
1、 算法的设计方案： （一）、总体方案设计：
（1）解非线性方程组。将给定的当作已知量代入题目给定的非线性方程组，求得与相对应 的数组t[i][j],u[i][j]。 （2）分片二次代数插值。通过分片二次代数插值运算，得到与数组t[11][21],u[11][21]] 对应的数组z[11][21]，得到二元函数z=。 （3）曲面拟合。利用x[i],y[j],z[11][21]建立二维函数表，再根据精度的要求选择适当k 值，并得到曲面拟合的系数矩阵C[r][s]。 （4）观察和的逼近效果。观察逼近效果只需要重复上面（1）和（2）的过程，得到与新的 插值节点对应的，再与对应的比较即可，这里求解可以直接使用（3）中的C[r][s]和k。
{ temp=0; for(l=k+1;l<=3;l++) {temp=temp+dF[k][l]*dx[l]/dF[k][k];} dx[k]=-F[k]/dF[k][k]-temp; } temp=0; for(l=0;l<=3;l++) /*求解矩阵范数，用无穷范数*/ { if(temp<fabs(dx[l])) temp=fabs(dx[l]); } fx=temp; temp=0; for(l=0;l<=3;l++) { if(temp<fabs(X[l])) temp=fabs(X[l]); } fX=temp; if(fabs(fx/fX)<Epsilon1) /*判断是否成立*/ { t[i][j]=X[0]; u[i][j]=X[1]; goto loop4;} else { for(l=0;l<=3;l++) {X[l]=X[l]+dx[l];} n=n+1; goto loop3;} } loop3:{if(n<M) /*判断是否超出规定迭代次数*/ goto loop1; else printf("迭代不成功\n"); goto loop4; } loop4:{continue;} } } } void fpeccz(double t[11][21],double u[11][21])/*分片二次代数插值子程序*/ { int s[11][21],r[11][21]; int i,j,i1,j1,m; double z0[6][6]={{-0.5,-0.34,0.14,0.94,2.06,3.5}, {-0.42,-0.5,-0.26,0.3,1.18,2.38}, {-0.18,-0.5,-0.5,-0.18,0.46,1.42},
对上面两式进行变形，得到如下两个线性方程组： ， 通过解上述两个线性方程组，则有： 3） 对于每一个， 。 4） 拟合需要达到的精度条件为： 。 其中对应着插值得到的数表中的值。 5） 让k逐步增加，每一次重复执行以上几步，直到 成立。此时的k值就是要求解最小的k。
2、 源程序：
#include<stdio.h> #include<iostream> #include<stdlib.h> #include<math.h> #include<float.h> #include<iomanip> #define Epsilon1 1e-12 /*解线性方程组时近似解向量的精度*/ #define M 200 /*解线性方程组时的最大迭代次数*/ #define N 10 /*求解迭代次数时假设的k的最大值，用于定义包含k的存储空间*/ void Newton(); /*牛顿法求解非线性方程组子程序*/ void fpeccz(); /*分片二次代数插值子程序*/ void qmnh(); /*曲面拟合子程序*/ void duibi(); /*对比 和p逼近效果的子程序*/ double x[11],y[21],t[11][21],u[11][21];/*定义全局变量*/ double z[11][21],C[10][10]; double kz; void Newton(double x[11],double y[21])/*牛顿法求解非线性方程组子程序*/ { double X[4],dx[4],F[4],dF[4][4],temp,m,fx,fX; int i,j,k,l,p,ik,n; for(i=0;i<=10;i++) { for(j=0;j<=20;j++) { X[0]=1; /*选取迭代初始向量，四个分别代表t，u，v，w*/ X[1]=1; X[2]=1; X[3]=1; n=0; loop1:{ F[0]=0.5*cos(X[0])+X[1]+X[2]+X[3]-x[i]-2.67; F[1]=X[0]+0.5*sin(X[1])+X[2]+X[3]-y[j]-1.07; F[2]=0.5*X[0]+X[1]+cos(X[2])+X[3]-x[i]-3.74; F[3]=X[0]+0.5*X[1]+X[2]+sin(X[3])-y[j]-0.79; /*求解F（x）*/