BUAA数值分析大作业三
北航数值分析大作业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 。
北航数值分析大作业三
一、题目：关于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
BUAAOS——Lab3实验报告
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)
北航数值分析报告大作业第三题(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。
北航数值分析大作业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。
中科院研究生院信息工程学院课件数值分析数值分析第三次作业及答案
中科院研究⽣院信息⼯程学院课件数值分析数值分析第三次作业及答案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)⼗。
北航《数值分析》习题
北航《数值分析》习题习题一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 。
北航数值分析大作业题目三
《数值分析》第三次大作业一、算法的设计方案： （一）、总体方案设计：（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 的效果。
数值分析大作业三、四、五、六、七
大作业 三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)<epflag=0;endx0=x1;endfprintf('方程的一个近似解为：%f\n',x0);寻找最大δ值的程序：cleareps=input('请输入搜索精度：');ep=input('请输入容许误差：');flag=1;k=0;x0=0;while flag==1sigma=k*eps;x0=sigma;k=k+1;m=0;flag1=1;while flag1==1 && m<=10^3x1=x0-fu(x0)/dfu(x0);if abs(x1-x0)<epflag1=0;endm=m+1;x0=x1;endif flag1==1||abs(x0)>=epflag=0;endend fprintf('最大的sigma 值为：%f\n',sigma);2．求下列方程的非零根5130.6651()ln 05130.665114000.0918x x f x x +⎛⎫=-= ⎪-⨯⎝⎭解：Matlab 程序为：（1）主程序clearclcformat longx0=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)>errorlimn=n+1;elsebreak ;endx0=x;enddisp(['迭代次数: n=',num2str(n)])disp(['所求非零根: 正根x1=',num2str(x),' 负根x2=',num2str(-x)])（2）子函数非线性函数ffunction y=f(x)y=log((513+*x)/*x))-x/(1400*;end（3）子函数非线性函数的一阶导数dffunction y=df()syms x1y=log((513+*x1)/*x1))-x1/(1400*;y=diff(y);end运行结果如下：迭代次数: n=5所求非零根: 正根x1= 负根x2=大作业 四试编写MATLAB 函数实现Newton 插值，要求能输出插值多项式. 对函数21()14f x x=+在区间[-5，5]上实现10次多项式插值.分析：（1）输出插值多项式。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
北京航空航天大学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四、讨论分析在编写程序过程中，前两步利用程序求解非线性方程组的过程中，可采用两种方式求解，一种是列消元法，另一种是牛顿迭代法求解非线性方程组。