下载文档原格式
数值分析第三次作业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 。
一、题目：关于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)。
BUAAOS——Lab3实验报告lab3实验报告思考题3.1 为什么我们在构造空闲进程链表时必须使⽤特定的插⼊的顺序？(顺序或者逆序)完成空闲链表的插⼊后，envs数组下标正好对应链表中的由前到后的顺序，因此调⽤空闲进程数组时优先调⽤下标最⼩的。
3.2 思考env.c/mkenvid 函数和envid2env 函数:请你谈谈对mkenvid函数中⽣成id 的运算的理解，为什么这么做？根据进程数组下标给每个进程⼀个独⼀⽆⼆的进程ID，同时还可以通过进程号获取数组偏移量找到进程块进⽽获得进程的全部信息。
为什么envid2env中需要判断e->env_id != envid 的情况？如果没有这步判断会发⽣什么情况？进程块也许会被替换，⽽ID只与⼀个进程块对应，如果没有这⼀步，可能会执⾏错误的进程块。
3.3 结合include/mmu.h中的地址空间布局，思考env_setup_vm 函数：我们在初始化新进程的地址空间时为什么不把整个地址空间的 pgdir都清零，⽽是复制内核的boot_pgdir作为⼀部分模板？(提⽰:mips 虚拟空间布局)因为每个进程都需要拥有内核的页表信息，这样才能够在陷⼊内核时正确执⾏，因此需要把内核拥有的虚拟地址对应的页表进⾏拷贝。
UTOP 和 ULIM 的含义分别是什么，在 UTOP 到 ULIM 的区域与其他⽤户区相⽐有什么最⼤的区别？UTOP是⽤户进程读写部分的最⾼地址，ULIM是⽤户进程的最⾼地址，UTOP到ULIM的区域为⽤户进程的进程块还有页表，不能被⽤户⾃⼰更改只能读取。
在 step4 中我们为什么要让pgdir[PDX(UVPT)]=env_cr3?(提⽰: 结合系统⾃映射机制)env_cr3是页⽬录的物理地址，pgdir[PDX(UVPT)]对应页⽬录的页⽬录项。
谈谈⾃⼰对进程中物理地址和虚拟地址的理解。
进程对应的地址都是虚拟地址，所有的物理地址只能从页表中查询所得到。
北航数值分析全部三次大作业第一次大作业是关于解线性方程组的数值方法。
我们被要求实现各种常用的线性方程组求解算法，例如高斯消元法、LU分解法和迭代法等。
我首先学习了这些算法的原理和实现方法，并借助Python编程语言编写了这些算法的代码。
在实验中，我们使用了不同规模和条件的线性方程组进行测试，并比较了不同算法的性能和精度。
通过这个作业，我深入了解了线性方程组求解的原理和方法，提高了我的编程和数值计算能力。
第二次大作业是关于数值积分的方法。
数值积分是数值分析中的重要内容，它可以用于计算曲线的长度、函数的面积以及求解微分方程等问题。
在这个作业中，我们需要实现不同的数值积分算法，例如矩形法、梯形法和辛普森法等。
我学习了这些算法的原理和实现方法，并使用Python编写了它们的代码。
在实验中，我们计算了不同函数的积分值，并对比了不同算法的精度和效率。
通过这个作业，我深入了解了数值积分的原理和方法，提高了我的编程和数学建模能力。
第三次大作业是关于常微分方程的数值解法。
常微分方程是数值分析中的核心内容之一，它可以用于描述众多物理、化学和生物现象。
在这个作业中，我们需要实现不同的常微分方程求解算法，例如欧拉法、龙格-库塔法和Adams法等。
我学习了这些算法的原理和实现方法，并使用Python编写了它们的代码。
在实验中，我们解决了一些具体的常微分方程问题，并比较了不同算法的精度和效率。
通过这个作业，我深入了解了常微分方程的原理和方法，提高了我的编程和问题求解能力。
总的来说，北航数值分析课程的三次大作业非常有挑战性，但也非常有意义。
通过这些作业，我在数值计算和编程方面得到了很大的提升，也更加深入地了解了数值分析的理论和方法。
虽然这些作业需要大量的时间和精力，但我相信这些努力将会对我未来的学习和工作产生积极的影响。
北航数值分析报告大作业第三题(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 SY1415215DIMENSIONX(11),Y(21),T(6),U(6),Z(6,6),UX(11,21),TY(11,21),FXY(11,21), C(6,6) DIMENSIONX1(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)) ENDDO ENDDO!*************分片二次插值得到z=f（x，y）***********DO I=1,11DO J=1,21CALL INTERPOLATION(UX(I,J),TY(I,J),FXY(I,J))ENDDO ENDDOWRITE (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(ETA).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)ENDIF ENDDODO 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(UX(N-1)-H/2) THENK=N-1ELSEDO I=3,N-2IF((U>X(I)-H/2).AND.(UY(M-1)-T/2) THENR=M-1 ELSEDO J=3,M-2IF((V>Y(J)-T/2).AND.(VN) P=N IF(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 ENDDO APX(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)*G ENDDOV(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=G2 G2=GENDDOD2=D2+G*GAPX(K)=APX(K)+X(I)*G*GDO J=1,M V(K,J)=V(K,J)+Z(I,J)*G ENDDOENDDODO 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) ENDDO U(J,1)=U(J,1)/D1ENDDOIF(Q>1)THEND2=0.0APY(2)=0.0DO I=1,MG=Y(I)-APY(1)D2=D2+G*G APY(2)=APY(2)+(Y(I))*G*G ENDDO APY(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=G2 G2=GENDDOD2=D2+G*GAPY(K)=APY(K)+Y(I)*G*G DO J=1,PU(J,K)=U(J,K)+V(J,I)*G ENDDOENDDODO 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) ENDDO ENDIFV(I,1)=-APY(I-1)*V(I-1,1)-BY(I-1)*V(I-2,1)ENDDO DO 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.0790.639 -0.18E+00 0.00 1.250 1.064 0.669 -0.52E+00 0.00 1.3001.050 0.699 -0.19E+00 0.00 1.350 1.037 0.729 -0.48E+00 0.001.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.2620.512 -0.11E+00 0.16 1.100 1.246 0.541 -0.74E+00 0.16 1.1501.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.2250.03E+000.24 0.650 1.482 0.275 0.64E+00 0.24 0.700 1.460 0.3010.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+00 0.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.3630.534 -0.25E+00 0.32 1.200 1.347 0.563 -0.29E+00 0.32 1.2501.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.1990.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.2230.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.4940.599 -0.82E+00 0.56 1.400 1.481 0.627 -0.26E+00 0.56 1.4501.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.1680.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+00 0.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.300 1.630 0.539 0.49E+00 0.72 1.350 1.616 0.5680.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.1510.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.4430.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-01 0.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+00 0.700 0.700 0.89E+01 0.29E+01 0.700 0.900 0.91E+01 0.11E+010.700 1.100 0.60E+00 0.97E+00 0.700 1.300 0.22E-01 0.06E-01 0.7001.500 -0.53E+00 -0.80E+00 0.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。
《数值分析》作业三院系：机械学院学号：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。
中科院研究⽣院信息⼯程学院课件数值分析数值分析第三次作业及答案6数值分析第三次作业及答案明当h T 0时，它收敛于原初值问题的准确解 y证：梯形公式为 y n ⼗ yn+—[f(X n ，y n )+f(X n^1，y n 』] h由 f (X,y) = —y= y n+ =y n +2( — y n — yn G=L n = 3 Y y n」訓 /乂⼚％l 2+h <12 +h ⼃丫2. (P202(6))写出⽤四阶经典的龙格⼀库塔⽅法求解下列初值问题的计算公式:y n + =y n + — (k 1 +2k 2 +2k 3 +k 4)=0.2214x n +1.2214y n +0.0214 6飞=3y n ⼼+x n )2)” k 2 =3(y n +0.1k 1〃(1+Xn +0.1) )L s =3(yn +0.1k2”(1+Xn +0.1)k 4 =3(yn +0.2k 3”(1+Xn +0.2) 0 2yn + =y n +〒(k 1 +2k 2 +2k 3 +k 4).1. ( P201 (4))⽤梯形⽅法解初值问题〔爲证明其近似解为y n 偌〕:并证⽤ . f 2-h 1 因 yoT " yn F ⼃.⽤上述梯形公式以步长h 经n 步计算到y n ,故有nh ：=x.X◎ T 茹Jf 2—h \n7 l 2+h ⼃1) ]y =x + y, 0 e x ￡1; ly(0) =1;2)l y \3%+x)，O *1； [y(0)=1.解：令h =0.2k 1 = f (X n , y n )= h k2=f (Xn+；；,yn+-k1)=Xn+- + 2 2 2 h k s = f (X n +；, y n +-k 2)=X n +- +y n +-k 2 =1.11(X n + y n )+0.11 2 2 2 2X n +y nh 1)4h h ??yn +；；k i =1.1(Xn +y n )+0.1 2 2 h . . h .................................. .2 ⼋ 2 J 'k 4 = f(X n +h,y n +hk 3)=X n + h + y n +hk 3 =1.222(X n +y n )+0.2223. (P202(7))证明对任意参数t，下列龙格库塔—公式是⼆阶的:r hy n 卄yn+^g+G)；* K i = f (X n, yj；K2 = f (X n +th, y n +thK i);[K3 =f(Xn+(1—t)h,yn+(1—t)hK i).证：由⼀元函数的泰勒展开有2 '''"y(X nG =y(X n) +hy'(X n)⼸[f x(X n,y(X n)) +f y(X n,y(X n))f(X n,y(X n))]中严h'2 3!⼜由⼆元函数的泰勒展开有y n41 =y n +；2(⼼+K3)=y n +；2[(f(X n,y n) + ￡%区『)也+f y(X n, y n)thf (X n, Y n)⼗。
北京航空航天大学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 相同。
将),(**j i y x 代入原方程组，求解响应),(**ij ij u t 进行分片双二次插值求得),(**j i y x f 。
),(**ji y x p 的计算则可以直接将),(**j i y x 代入所求p(x,y)。
二、 源程序********* 第三次数值分析大作业Q9************ integer::i, j, K1, L1, n, mdimension X(31), Y(21), T(6), U(6), Z(6, 6), UX(11, 21), TY(11, 21), FXY(11, 21), C(6, 6) dimension z1(31, 21), z2(31, 21), z3(31, 21), z4(31, 21), z5(31, 21) 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, Error, X1, Y1, FXY1, PXY1, UX1, TY1, z1, z2, z3, z4, z5OPEN(1, FILE = '数值分析Q9程序结果.TXT')do i = 1, 31X(i) = 1 + 0.008 * (i - 1)end dodo j = 1, 21Y(j) = 1 + 0.008 * (j - 1) !生成矩阵end do!*****求解非线性方程组，得到数表(X(i), Y(j), f1, f2, f3, f4, f5)*******do i = 1, 31do j = 1, 21call Newton(X(i), Y(j), z1(i, j), z2(i, j), z3(i, j), z4(i, j), z5(i, j)) !用离散牛顿法求解非线性方程组。
if (i == 27.and.j == 21)then!write(*, *)z1(i, j), z2(i, j), z3(i, j), z4(i, j), z5(i, j)!pauseend ifend doend dowrite(*, "('数表( X(i),Y(j),f1,f2,f3,f4,f5 ):')")write(*, "(3X,'X',7X,'Y',10X,'f1(x,y)',13x,'f2(x,y)',13x,'f3(x,y)',13x,'f4(x,y)',13x,'f5(x,y)')")do i = 1, 31do j = 1, 21write(1, '(1X,F5.3,2X,F5.3,2X,5E20.13)') X(i), Y(j), z1(i, j), z2(i, j), z3(i, j), z4(i, j), z5(i, j) write(*, '(3X,I2,1X,I2,3X,F5.3,2X,F5.3,2X,5E20.13)')i, j, X(i), Y(j), z1(i, j), z2(i, j), z3(i, j), z4(i, j), z5(i, j)end doend do!***********最小二乘拟合得到P（x，y）**************N = 31M = 21WRITE(1, '(" ","K和σ分别为：")')WRITE(*, '(" ","选择过程的K1,L1和σ分别为：")')DO K1 = 1, 8DO L1 = 1, 8CALL surface_fit(X, Y, Z4, C, M, N, K1, L1, Error) !surface_fit(X, Y, FXY, C, N, M, K, K, E) K是最高次项，E是误差，C是系数矩阵WRITE(*, '(I3,2X,I3,2X,E20.13)') K1 - 1, L1 - 1, Error !subroutinesurface_fit(X, Y, Z, A, N, M, P, Q, DT1)IF(Error<0.10e-4 . or . (K1 == 5. and .L1 == 4))thenWRITE(*, '(" 目前精度符合给定要求")')goto 30end ifpauseENDDOENDDO30 continueWRITE(1, '(" ")')WRITE(*, '("Q9的系数矩阵Crs（按行）为：")')DO I = 1, K1DO J = 1, L1WRITE(*, '(E20.13,2X)') C(I, J)ENDDOWRITE(1, "('')")WRITE(*, "('')")ENDDOpauseCLOSE(1)END!***********用离散牛顿法求解非线性方程组****************subroutine Newton(X1, Y1, z1, z2, z3, z4, z5) !返回z1, z2, z3, z4, z5PARAMETER(N = 5)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, E2, z1, z2, z3, z4, z5F1(z1, z2, z3, z4, z5) = exp(-z1) + exp(-2 * z2) + z3 - 2 * z4 + X1 * z5 - 5.3F2(z1, z2, z3, z4, z5) = exp(-2 * z1) + exp(-z2) - 2 * z3 + Y1 * z4 - z5 + 25.6F3(z1, z2, z3, z4, z5) = X1 * z1 + 3 * z2 + exp(-z3) - 3 * z5 + 37.8F4(z1, z2, z3, z4, z5) = 2 * z1 + Y1 * z2 + z3 - exp(-z4) + 2 * exp(-2 * z5) - 31.3F5(z1, z2, z3, z4, z5) = z1 - 2 * z2 - 3 * X1 * z3 + exp(-2 * z4) + 3 * exp(-z5) + 42.1EPS = 10E-15 !迭代误差限制M = 100 !最大迭代次数X = 2.0 !x为初始迭代向量do K = 1, MH = 2 !求差商的步长不可以取做0.1!计算Y = F（x）Y(1) = F1(X(1), X(2), X(3), X(4), X(5))Y(2) = F2(X(1), X(2), X(3), X(4), X(5))Y(3) = F3(X(1), X(2), X(3), X(4), X(5))Y(4) = F4(X(1), X(2), X(3), X(4), X(5))Y(5) = F5(X(1), X(2), X(3), X(4), X(5))!计算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)end if!write(*, *)i, j, e(jj)end do!pauseif (I == 1) thenJA(I, J) = (F1(E(1), E(2), E(3), E(4), E(5)) - F1(X(1), X(2), X(3), X(4), X(5))) / H(J) elseif(I == 2) thenJA(I, J) = (F2(E(1), E(2), E(3), E(4), E(5)) - F2(X(1), X(2), X(3), X(4), X(5))) / H(J) elseif(I == 3) thenJA(I, J) = (F3(E(1), E(2), E(3), E(4), E(5)) - F3(X(1), X(2), X(3), X(4), X(5))) / H(J) elseif(I == 4) thenJA(I, J) = (F4(E(1), E(2), E(3), E(4), E(5)) - F4(X(1), X(2), X(3), X(4), X(5))) / H(J) elseif(I == 5) thenJA(I, J) = (F5(E(1), E(2), E(3), E(4), E(5)) - F5(X(1), X(2), X(3), X(4), X(5))) / H(J) end if!write(*, *)i, j, JA(I, J)!pauseend doend do!求解线性方程组CALL LU(JA, XK, -Y, N) !GAUSS(A, X, B, N) 返回X AX = B!判断精度CALL NORM(XK, N, E1) !NORM(X, N, A)求向量X的无穷范数，即求向量X中最大值A 返回值为最大值ACALL NORM(X, N, E2)if (E1 / E2 <= EPS) thenz1 = X(1)z2 = X(2)z3 = X(3)z4 = X(4)z5 = X(5)exitelsedo I = 1, NX(I) = X(I) + XK(I)!write(*, *)k, i, X(i), xk(i), E1 / E2end doend ifend do!write(*, *)k, '结束'!pauseRETURNEND!**********利用列主元高斯消去法求解线性方程组* ********subroutine GAUSS(A, X, B, N) !GAUSS(A, X, B, N) 返回X AX = Bdimension 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)end doTB(K) = B(K)B(K) = B(TL)B(TL) = TB(K)end ifend dodo 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)end doB(I) = B(I) - M(I, K) * B(K)end doend do!回代过程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)end doX(K) = (B(K) - S) / A(K, K)end doRETURNENDsubroutine LU(aa, x, b, n) !GAUSS(A, X, B, N) 返回X AX = B implicit noneinteger(4)::n, sgnreal(8)::aa(n, n), b(n), x(n)integer(4)::r, i, k, imaxreal(8), allocatable::a(:, : )real(8)::tallocate(a(n, n))a = aa;sgn = 1!分解aa = LU，L和U存放在a中do r = 1, n!列主元过程do i = r, ndo k = 1, r - 1a(i, r) = a(i, r) - a(i, k) * a(k, r)end doend doimax = rdo i = r + 1, nif (abs(a(i, r)) > abs(a(imax, r))) imax = i end doif (abs(a(imax, r)) == 0) thensgn = 0;EXITend if!交换第r行和第imax行if (imax /= r) thendo i = 1, nt = a(r, i);a(r, i) = a(imax, i);a(imax, i) = t end dot = b(r);b(r) = b(imax);b(imax) = tend if!计算l的第r列和u的第r行do i = r + 1, na(i, r) = a(i, r) / a(r, r)do k = 1, r - 1a(r, i) = a(r, i) - a(r, k) * a(k, i)end doend doend do!sgn = 1说明分解成功，可继续求解if (sgn == 1) then!求解Ly = b，y放入x中do i = 1, nx(i) = b(i)do k = 1, i - 1x(i) = x(i) - a(i, k) * x(k)end doend do!求解Ux = ydo i = n, 1, -1do k = n, i + 1, -1x(i) = x(i) - a(i, k) * x(k)end dox(i) = x(i) / a(i, i)end doend ifend subroutine LU!***********求向量的无穷范数************subroutine NORM(X, N, A) !NORM(X, N, A)求向量X的无穷范数，即求向量X中最大值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))end ifend doRETURNEND!**************分片二次代数插值**************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 = Iend ifend doend ifif (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 = Jend ifend doend ifI = 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)end doend doRETURNEND!*******************最小二乘法拟合子函数************** subroutine surface_fit(X, Y, Z, A, M, N, P, Q, DT1)INTEGER P, Qdimension X(N), Y(M), Z(N, M), A(P, Q)dimension APX(30), APY(30), BX(30), BY(30), U(30, 20), V(30, M) dimension T(20), T1(20), T2(20)real(8) X, Y, Z, A, DT1do I = 1, Ndo J = 1, M!write(*, *)X(I), Y(J), Z(I, J)end doend dodo I = 1, Pdo J = 1, QA(I, J) = 0.0end doend doif (P > N) P = Nif (P > 8) P = 8if (Q > M) Q = Mif (Q > 8) Q = 8XX = 0YY = 0D1 = NAPX(1) = 0.0do I = 1, NAPX(1) = APX(1) + X(I)end doAPX(1) = APX(1) / D1do J = 1, MV(1, J) = 0.0do I = 1, NV(1, J) = V(1, J) + Z(I, J)end doV(1, J) = V(1, J) / D1end doif (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 end doAPX(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) * Gend doV(2, J) = V(2, J) / D2end doD1 = D2end ifdo K = 3, PD2 = 0.0APX(K) = 0.0do J = 1, MV(K, J) = 0.0end dodo 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 = Gend doD2 = D2 + G * GAPX(K) = APX(K) + X(I) * G * G do J = 1, MV(K, J) = V(K, J) + Z(I, J) * G end doend dodo J = 1, MV(K, J) = V(K, J) / D2end doAPX(K) = APX(K) / D2BX(K) = D2 / D1D1 = D2end doD1 = MAPY(1) = 0.0do I = 1, MAPY(1) = APY(1) + Y(I)end doAPY(1) = APY(1) / D1do J = 1, PU(J, 1) = 0.0do I = 1, MU(J, 1) = U(J, 1) + V(J, I)end doU(J, 1) = U(J, 1) / D1end doif (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 end doAPY(2) = APY(2) / D2BY(2) = D2 / D1do J = 1, PU(J, 2) = 0.0G = Y(I) - APY(1)U(J, 2) = U(J, 2) + V(J, I) * Gend doU(J, 2) = U(J, 2) / D2end doD1 = D2end ifdo K = 3, QD2 = 0.0APY(K) = 0.0do J = 1, PU(J, K) = 0.0end dodo 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 = Gend doD2 = D2 + G * GAPY(K) = APY(K) + Y(I) * G * Gdo J = 1, PU(J, K) = U(J, K) + V(J, I) * Gend doend dodo J = 1, PU(J, K) = U(J, K) / D2end doAPY(K) = APY(K) / D2BY(K) = D2 / D1D1 = D2end doV(1, 1) = 1.0V(2, 1) = -APY(1)V(2, 2) = 1.0do I = 1, Pdo J = 1, Qend doend dodo 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) end doend ifV(I, 1) = -APY(I - 1) * V(I - 1, 1) - BY(I - 1) * V(I - 2, 1)end dodo 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)end doend ifT(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)end doend ifdo J = 1, Qdo K = I, 1, -1do L = J, 1, -1A(K, L) = A(K, L) + U(I, J) * T(K) * V(J, L) end doend doend doend doDT1 = 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)end doG = G * X2DD = DD + GX2 = X2 * X1end doDT = DD - Z(I, J)DT1 = DT1 + DT * DTend doend doRETURNEND三、结果所有Z1值输出：当t=1时,k=1,l=1时，误差为8.837746002099e-001 当t=1时,k=1,l=2时，误差为6.977120426215e-001 当t=1时,k=1,l=3时，误差为6.975657365251e-001 当t=1时,k=1,l=4时，误差为6.975471855564e-001 当t=1时,k=1,l=5时，误差为6.975459514363e-001 当t=1时,k=1,l=6时，误差为6.975460467607e-001 当t=1时,k=1,l=7时，误差为6.975458615969e-001 当t=1时,k=1,l=8时，误差为6.975458438559e-001 当t=1时,k=2,l=1时，误差为1.898010561451e-001 当t=1时,k=2,l=2时，误差为3.280174917380e-003 当t=1时,k=2,l=3时，误差为3.046761146841e-003 当t=1时,k=2,l=4时，误差为3.019691796469e-003 当t=1时,k=2,l=5时，误差为3.017656567704e-003 当t=1时,k=2,l=6时，误差为3.017822313207e-003 当t=1时,k=2,l=7时，误差为3.017501224160e-003 当t=1时,k=2,l=8时，误差为3.017470387714e-003 当t=1时,k=3,l=1时，误差为1.876627604962e-001当t=1时,k=3,l=3时，误差为6.274729155708e-004 当t=1时,k=3,l=4时，误差为5.958207057060e-004 当t=1时,k=3,l=5时，误差为5.932790508382e-004 当t=1时,k=3,l=6时，误差为5.934981393173e-004 当t=1时,k=3,l=7时，误差为5.930751842016e-004 当t=1时,k=3,l=8时，误差为5.930344418470e-004 当t=1时,k=4,l=1时，误差为1.872704724551e-001 当t=1时,k=4,l=2时，误差为4.391840464527e-004 当t=1时,k=4,l=3时，误差为1.553040890731e-004 当t=1时,k=4,l=4时，误差为1.217532961144e-004 当t=1时,k=4,l=5时，误差为1.189542633273e-004 当t=1时,k=4,l=6时，误差为1.192056698476e-004 当t=1时,k=4,l=7时，误差为1.187214970005e-004 当t=1时,k=4,l=8时，误差为1.186747956604e-004 当t=1时,k=5,l=1时，误差为1.871946647750e-001 当t=1时,k=5,l=2时，误差为3.472455970062e-004 当t=1时,k=5,l=3时，误差为5.962555041392e-005 当t=1时,k=5,l=4时，误差为2.539055014905e-005 当t=1时,k=5,l=5时，误差为2.247865931059e-005 当t=1时,k=5,l=6时，误差为2.274752542131e-005 当t=1时,k=5,l=7时，误差为2.223142517293e-005 当t=1时,k=5,l=8时，误差为2.218158914911e-005 当t=1时,k=6,l=1时，误差为1.871815109880e-001 当t=1时,k=6,l=2时，误差为3.305967713992e-004 当t=1时,k=6,l=3时，误差为4.197780459490e-005 当t=1时,k=6,l=4时，误差为7.517279189052e-006 当t=1时,k=6,l=5时，误差为4.560806407266e-006 当t=1时,k=6,l=6时，误差为4.854134849636e-006 当t=1时,k=6,l=7时，误差为4.309636811778e-006 当t=1时,k=6,l=8时，误差为4.255536604730e-006 当t=1时,k=7,l=1时，误差为1.871811559512e-001 当t=1时,k=7,l=2时，误差为3.301355384546e-004 当t=1时,k=7,l=3时，误差为4.148264324908e-005 当t=1时,k=7,l=4时，误差为7.013642621233e-006 当t=1时,k=7,l=5时，误差为4.0558********e-006 当t=1时,k=7,l=6时，误差为4.332957246970e-006 当t=1时,k=7,l=7时，误差为3.803282499915e-006 当t=1时,k=7,l=8时，误差为3.750510828545e-006 当t=1时,k=8,l=1时，误差为1.871814689902e-001 当t=1时,k=8,l=2时，误差为3.305647770605e-004 当t=1时,k=8,l=3时，误差为4.195636990633e-005当t=1时,k=8,l=5时，误差为4.544871301732e-006当t=1时,k=8,l=6时，误差为4.823190957063e-006当t=1时,k=8,l=7时，误差为4.294113913876e-006当t=1时,k=8,l=8时，误差为4.240361624684e-0061.达到精度要求时的k t,l t,σt值以及p t(x,y)中的系数c rs(r=0,1,⋯,k t;s=0,1,⋯,l t;t=1,2,3,4,5)当t=1时,k=5,l=3满足精度要求，误差为5.962555041392e-005此时C_rs系数为：1.441333590621e+006 -4.116031068916e+006 3.916977150672e+006 -1.242220462015e+006-6.590516181073e+006 1.882151198097e+007 -1.791206005171e+007 5.680829583415e+0061.204357289289e+007 -3.439593845184e+007 3.273522974877e+007 -1.038246720132e+007-1.0994********e+007 3.140123711855e+007 -2.988636994850e+007 9.479358835974e+0065.014097270085e+006 -1.432119726471e+007 1.363093714658e+007 -4.323683876383e+006-9.138966612586e+005 2.610377236358e+006 -2.484685814172e+006 7.881771039422e+005 当t=2时,k=6,l=4满足精度要求，误差为9.0353********e-005此时C_rs系数为：-5.481891067578e+008 2.074011382819e+009 -2.941534940226e+0091.853*********e+009 -4.378840934629e+0082.998232943450e+009 -1.134375350415e+010 1.608907750759e+010 -1.013875217844e+010 2.395200867865e+009-6.826947365615e+009 2.583028501991e+010 -3.663663812558e+010 2.308776682663e+010 -5.454479189947e+0098.283724891885e+009 -3.134296670524e+010 4.445689888177e+010 -2.801683420792e+010 6.619192714823e+009-5.649221682013e+009 2.137545840994e+010 -3.0319********e+010 1.910829477601e+010 -4.514640958556e+0092.0530********e+009 -7.768443759281e+009 1.101946026381e+010 -6.944955070352e+009 1.640917711621e+009-3.106243267156e+008 1.175408955519e+009 -1.667362212117e+009 1.050884030462e+009 -2.483070415140e+008当t=3时,k=5,l=5满足精度要求，误差为9.977362178063e-005此时C_rs系数为：3.600089486845e+008 -1.699533575604e+009 3.208447729356e+009 -3.027*********e+009 1.428353462225e+009 -2.694766123536e+008-1.640687195429e+009 7.745516701144e+009 -1.462259337291e+010 1.379958686345e+010 -6.510057311810e+009 1.228234336256e+0092.987753699081e+009 -1.410516149964e+010 2.662940013115e+010 -2.513119188736e+010 1.185612079955e+010 -2.236921263616e+009-2.717575309324e+009 1.282992245959e+010 -2.422240189463e+010 2.286017284628e+010 -1.078500769867e+010 2.034889642746e+0091.234633202632e+009 -5.828947079425e+009 1.100509309496e+010 -1.038645262295e+010 4.900272955242e+009 -9.245995848846e+008-2.241332456279e+008 1.058202394055e+009 -1.997944415821e+009 1.885683504372e+009 -8.896815924505e+008 1.678733636367e+008当t=4时,k=4,l=3满足精度要求，误差为7.394179889802e-005此时C_rs系数为：8.152974907438e+004 -2.345663684336e+005 2.250345834835e+005 -7.200899226698e+004-2.993600430891e+005 8.612568081229e+005 -8.263127337022e+005 2.644409818494e+0054.117428501866e+005 -1.184513770565e+006 1.136583298861e+006 -3.637833993098e+005-2.513640576849e+005 7.232151659717e+005 -6.940789498845e+005 2.221891520521e+0055.748650032747e+004 -1.654317203876e+005 1.587997329300e+005 -5.084378431839e+004当t=5时,k=6,l=4满足精度要求，误差为6.973923411327e-005此时C_rs系数为：-4.813754261040e+008 1.821248788357e+009 -2.583068339065e+0091.627720921048e+009 -3.845267616752e+0082.632819230066e+009 -9.961308846347e+009 1.412843040695e+010 -8.903274694451e+009 2.103338621105e+009-5.994915737425e+009 2.268240438000e+010 -3.217200399140e+010 2.027*********e+010 -4.789821619187e+0097.274143893313e+009 -2.752321367855e+010 3.903914731265e+010 -2.460264378714e+010 5.812579350258e+009-4.960705502680e+009 1.877035952038e+010 -2.662481556699e+010 1.677961216602e+010 -3.964459600296e+0091.802796103644e+009 -6.821632608967e+009 9.676446317171e+009 -6.098539035943e+009 1.440932517943e+009-2.727624224446e+008 1.032142200109e+009 -1.464135924185e+009 9.227981346741e+008 -2.180424913973e+008四、讨论分析在编写程序过程中，前两步利用程序求解非线性方程组的过程中，可采用两种方式求解，一种是列消元法，另一种是牛顿迭代法求解非线性方程组。
