北航数值分析全部三次大作业

合集下载

北航数值分析大作业第一题幂法与反幂法

数值分析（B ）大作业（一）姓名：学号：电话：1、算法设计：①求1λ、501λ和s λ的值：s λ：s λ表示矩阵的按模最小特征值，为求得s λ直接对待求矩阵A 应用反幂法即可。

1λ、501λ：若矩阵A 的特征值满足关系 1n λλ<<且1n λλ≠，要求1λ、及501λ时，可按如下方法求解： a ．对矩阵A 用幂法，求得按模最大的特征值1m λ。

b ．按平移量1m λ对矩阵A 进行原点平移得矩阵1m BA I λ=+，对矩阵B 用反幂法求得B 的按模最小特征值2m λ。

c ． 321m m m λλλ=-则：113min(,)m m λλλ=，13max(,)n m m λλλ=即为所求。

②求和A 的与数5011140k k λλμλ-=+最接近的特征值ik λ（k=0，1，…39）：求矩阵A 的特征值中与P 最接近的特征值的大小，采用原点平移的方法：先求矩阵 B=A-PI 对应的按模最小特征值k β，则k β+P 即为矩阵A 与P 最接近的特征值。

在本次计算实习中则是先求平移矩阵k B A I μ=-，对该矩阵应用反幂法求得s λ，则与k μ最接近的A 的特征值为：s P λ+重复以上过程39次即可求得ik λ（k=0，1，…39）的值。

③求A 的（谱范数）条件数2cond()A 和行列式det A ：在（1）中用反幂法求矩阵A 的按模最小特征值时，要用到Doolittle 分解方法，在Doolittle 分解完成后得到的两个矩阵分别为L 和U ，则A 的行列式可由U 阵求出，即：det(A)=det(U)。

求得det(A)不为0，因此A 为非奇异的实对称矩阵，则： max 2()scond A λλ=，max λ和s λ分别为模最大特征值与模最小特征值。

2、程序源代码：#include "Stdio.h"#include "Conio.h"#include "math.h"//****************************************************************************// // 在存储带状矩阵时，下面的几个量在程序中反复用到，为方便编程故把它们定义成宏.// // M ：转换后的矩阵的行数，M=R+S+1。

北航数值分析大作业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 。

数值分析大作业三四五六七

数值分析大作业三四五六七数值分析大作业三四五六七Document number【SA80SAB-SAA9SYT-SAATC-SA6UT-SA18】大作业三1. 给定初值0x 及容许误差，编制牛顿法解方程f (x )=0的通用程序. 解：Matlab 程序如下：函数m 文件：fu.mfunction Fu=fu(x)Fu=x^3/3-x;end函数m 文件：dfu.mfunction Fu=dfu(x)Fu=x^2-1;end用Newton 法求根的通用程序Newton.mclear;x0=input('请输入初值x0:');ep=input('请输入容许误差:');flag=1;while flag==1x1=x0-fu(x0)/dfu(x0);if abs(x1-x0)<ep< p="">flag=0;endx0=x1;endfprintf('方程的一个近似解为：%f\n',x0);寻找最大δ值的程序：Find.mcleareps=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)endm=m+1;x0=x1;endif flag1==1||abs(x0)>=epflag=0;endendfprintf('最大的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<n< p="">x=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+0.6651*x)/(513-0.6651*x))-x/(1400*0.0918);end（3）子函数非线性函数的一阶导数dffunction y=df()syms x1y=log((513+0.6651*x1)/(513-0.6651*x1))-x1/(1400*0.0918);y=diff(y);end运行结果如下：迭代次数: n=5所求非零根: 正根x1=767.3861 负根x2=-767.3861大作业四试编写MATLAB 函数实现Newton 插值，要求能输出插值多项式. 对函数21()14f x x=+在区间[-5，5]上实现10次多项式插值.分析：（1）输出插值多项式。

数值分析第三次作业

26.解：(1).J 法：J ∴法收敛.GS 法：()11102210221101110122100210G B D L U ----⎡⎤⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=-=-=--⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦⎣⎦022023002-⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦()()()212322det 0232020,2,21G G I B B λλλλλλλλλρ--=-=--∴===∴=GS ∴法不收敛.()2.J 法：()()131231*********()12202101121101212012125det 412125550,,,1222J J J B D L U I B i i B λλλλλλλλλρ---⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥=+=--=--⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦--==+--∴===-∴=J ∴法不收敛.GS 法：()()()1312310220221101101122022022det 11002201J J J B D L U I B B λλλλλλλλρ---⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥=+=--=--⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥----⎣⎦⎣⎦⎣⎦--===⇒===∴=()1120111200011220212120021120014120G B D L U ----⎡⎤⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=-=-=--⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥--⎣⎦⎣⎦⎣⎦⎣⎦01212012120012-⎡⎤⎢⎥=--⎢⎥⎢⎥-⎣⎦()()()21231212det 01212120120,12,121G G I B B λλλλλλλλλρ--=+=+=+∴===-=GS ∴法收敛27.解：()1010911102,702106A b -⎛⎫⎛⎫ ⎪ ⎪=--= ⎪ ⎪⎪ ⎪-⎝⎭⎝⎭A 为严格对角占优阵，J ∴法和GS 法均收敛.J 法的分量形式为：()()()11111,1,2,,i nk k k ii ij j ij j j j i ii x b a x a x i n a -+==+⎛⎫=--= ⎪⎝⎭∑∑J ∴法的迭代格式为：(1)()12(1)()()213(1)()321(9)101(72)101(62)10k k k k k k k x x x x x x x +++⎧=+⎪⎪⎪=++⎨⎪⎪=+⎪⎩取初值(0)0x =，J 法的数值结果是：迭代次数k()1k x ()2k x ()3kx 1 0.900000 0.700000 0.600000 2 0.970000 0.910000 0.740000 3 0.991000 0.945000 0.782000 4 0.994000 0.955500 0.789000 50.9955500.9572500.791100GS 法的分量形式为：()()()111111,1,2,,i nk k k ii ij j ij j j j i ii x b a x a x i n a -++==+⎛⎫=--= ⎪⎝⎭∑∑GS ∴法的迭代格式为：(1)()12(1)(1)()213(1)(1)321(9)101(72)101(62)10k k k k k k k x x x x x x x +++++⎧=+⎪⎪⎪=++⎨⎪⎪=+⎪⎩取初值(0)0x =，GS 法的数值结果是：迭代次数k ()1k x()2k x()3kx10.900000 0.790000 0.758000 2 0.979000 0.949500 0.789900 3 0.994950 0.957475 0.791495 4 0.9957475 0.9578738 0.7915748 50.9957874 0.95789370.7915787()123210,99,950A∆=∆=∆=∴对称正定，()1110000101001011100102010105020110000105J B D L U -⎛⎫⎛⎫ ⎪ ⎪⎛⎫ ⎪⎪ ⎪⎪ ⎪=+=-= ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎝⎭ ⎪ ⎪ ⎪⎪⎝⎭⎝⎭212,310101111det()()00,10520201051()20J J I B B λλλλλλλλρ-∴-=--=-=⇒==±-∴=∴SOR 法的最优松弛因子为：[]2221.01282111/2011()opt J B ωρ==≈+-+-()10.01282opt opt L ρωω=-=对应的渐近收敛率为：R=-ln ()() 4.35654opt opt R L L ωρω=SOR 法的分量形式为：()()()()()111111,1,2,,i nk k k k iii ij j ijjj j i ii x x b a x a xi n a ωω-++==+⎛⎫=-+--= ⎪⎝⎭∑∑∴SOR 法（ω取最佳松弛因子）的迭代格式为：(1)()()112(1)()(1)()2213(1)()(1)3321.012820.01282(9)101.012820.01282(72)101.012820.01282(62)10k k k k k k k k k k x x x x x x x x x x +++++⎧=-++⎪⎪⎪=-+++⎨⎪⎪=-++⎪⎩取初值(0)0x =，SOR 法的数值结果是：迭代次数k ()1k x()2k x()3kx10.911538 0.801296 0.770006 2 0.981009 0.954035 0.791074 3 0.995588 0.957822 0.791571 4 0.995785 0.957894 0.791579 50.9957890.9578950.79157928.一定收敛.证明：对于11122122a a A a a ⎛⎫=⎪⎝⎭，A 对称正定，()212211122122111221201,2,,det()0,iia i a a A a a a a a a a ∴===-对于J 法：121111121121222221000()01000J a aa a B D L U a a a a -⎛⎫⎛⎫-⎪ ⎪-⎛⎫⎪ ⎪=+== ⎪ ⎪ ⎪-⎝⎭- ⎪ ⎪⎝⎭⎝⎭122211212121,2121122112222det()0J a a a a I B a a a a a a λλλλλ-==-=⇒=22121122121122()1J a a a a B a a ρ∴=∴J 法收敛. GS 法：12221112121221122112212112222121122121122det()00,0()1G G a a a a I B a a a a a a a a a a a B a a λλλλλλλρ⎛⎫-==-=⇒==⎪⎝⎭-∴=∴GS 法收敛.∴对于系数矩阵对称正定的2阶线性方程组，J 法和GS 法一定收敛. 30.证明：由线性代数知识知：∃可逆矩阵使121s J J p B P J J -⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎣⎦其中,i in n iJ R ⨯∈对应于特征值()121,2,,i s i s n n n nλ=++=()0B ρ=∴B 的所有特征值为0，120101,1,2,,10i i r r i s J R i s r r r n⨯⎛⎫⎪⎪⎪∴=∈=++= ⎪ ⎪ ⎪⎝⎭1i r =时，11i J R ⨯∈1i r 时， 0,i i r r i i J J R ⨯=∈12kkkk s J J J J ⎡⎤⎢⎥⎢⎥∴=⎢⎥⎢⎥⎢⎥⎣⎦最多迭代到第n 次，即k=n 时，10,0k k k J B PJ P -=== 设x *是Ax b =的精确解，误差向量()()k k e x x *=-()()()()()()110k k k k ke x x B x x B e B e--**=-=-===所以最多迭代到第n 次时，()()()00,k k k e B e x x *===所以结论成立31．解：（1）根据迭代公式(1)()()()k k k x x Ax b α+=--，有： (1)()()k k xI A x b αα+=-+ ∴迭代矩阵13212B I A ααααα--⎛⎫=-= ⎪--⎝⎭ 12132det()(1)(14)0121,14I B λααλλαλααλαλαλα-+∴-==-+-+=-+∴=-=-当{}()max 1,141B ραα=--时，迭代收敛111110121411141ααααα⎧---⎧⎪⇒⇒⎨⎨---⎪⎩⎩012α∴时，此迭代方法收敛{}()()m a x 1,141,00.441,0.40.5B B ρααααραα=---⎧∴=⎨-≤⎩ 0.4α∴=时，()Bρ最小，迭代收敛最快()12,n λλ为A 的特征值，11,1n αλαλ∴--为I A α-的特征值{}1()m a x 1,1n I A ρααλαλ∴-=--必要性：迭代收敛()1I A ρα⇒-111110211nαλαλαλ-⎧-⎪∴⇒⎨-⎪⎩所以必要性成立充分性：()1111102022,1,2,11,1,2,()max 11i i ini i i ni nI A αλαλλλαλρααλ--=∴≤=∴-=∴-=-所以此迭代法收敛，充分性成立（3） 1102αλ-时，111121,0()21,2n n in I A αλαλλρααλαλλλ-⎧-≤⎪+⎪-=⎨⎪-⎪+⎩根据图像，12nαλλ=+时，()I A ρα-最小33.解：()()()()()()()()()()()()()()()1()1121()111211111(1)()11111,k k k k k k x D L Ux D L b x L D U x L b x D U L D L Ux D U L D L b D U bC D U L D L U g D U L D L b D U b+--++-------+-----⎧=-+-⎪⎨⎪=--⎩⇒=--+--+-∴=--=--+-分析收敛性：()()()()()1111L D L D L DD LI D D L----=--+-=-+-⎡⎤⎣⎦()()()1111D D L I D L D D D L LD ----⎡⎤=---=-⎣⎦()()111C D U D D L LD U---∴=--()()()()()11112D L D D U L D U D L I DL D U L D U D L U-------=----=--=()()111I C I D L D D U LD U ---⎡⎤∴-=---⎣⎦()()()()111I D L D D U D L D D U A ---⎡⎤⎡⎤=------⎣⎦⎣⎦()()()()1111I I D L D D U A D L D D U A ----⎡⎤⎡⎤=-+--=--⎣⎦⎣⎦ 令()()1M D L D DU -⎡⎤=--⎣⎦1C I M A -∴=-因为A 对称正定，所以D 也正定令 1111222,()D D D W D D U ----==-TM W W ∴=()()11111112222TWCW D D L LD D D L LD -----⎡⎤⎡⎤=--⎢⎥⎢⎥⎣⎦⎣⎦ 令()11122P DD L LD--=-1T W C W P P -∴=所以C 与T PP 相似，其特征值均为非负实数1111()T I WCW W I C W WM AW W AW -------=-==所以 1I WCW --为对称正定矩阵，其特征值()110WCW λ--C ⇒的特征值()C λ满足()01C λ≤，故该迭代法收敛35.解：1112112111122212Ax Bx b x A Bx A b Bx Ax b x A Bx A b ----⎧+==-+⎧⎪⇒⎨⎨+==-+⎪⎩⎩∴J 法的迭代公式为：(1)()111111(1)()1222110000k k k k J x x A b A B A Bx x A b A B C A B +---+---⎛⎫⎛⎫⎛⎫⎛⎫-∴=+ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪-⎝⎭⎝⎭⎝⎭⎝⎭⎛⎫-∴= ⎪-⎝⎭若λ为矩阵1A B --的特征值，对应的特征向量为11111,0n x R x A Bx x λ-∈≠∴-= 11111111111111111111J J x x x A B x C x x x A B x x x x A B x C x x x A B x λλλλλλ----⎡⎤-⎡⎤⎡⎤⎡⎤===⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦⎣⎦-⎡⎤⎡⎤⎡⎤⎡⎤===-⎢⎥⎢⎥⎢⎥⎢⎥---⎣⎦⎣⎦⎣⎦⎣⎦∴若n 阶矩阵1A B --有特征值12,,,n λλλ，则2n 阶矩阵J C 有特征值12,,,nλλλ±±±38.（1）解：因为系数矩阵A 对称正定，所以可以运用共轭梯度法（CG ）解此方程组取()()00,0Tx =,()()()()000r p b Ax 0,1T∴==-=-，()()()()()()00r ,r 12p ,Ap α==()()()0,-T10001x xp2α⎛⎫∴=+= ⎪⎝⎭，()()(),0T10003r rAp2α⎛⎫=-=- ⎪⎝⎭,()()()()()(),11r ,r 94r r β==()()(),-T110039p r p 24β⎛⎫=+= ⎪⎝⎭()()()()()()11111r ,r 12p ,Ap α==,()()()()1,-2T2111x x p α=+=()2x 即为所求方程的精确解。

北航数值分析大作业三

一、题目：关于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)。

北航数值分析报告大作业第三题(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。

BUAA数值分析大作业三

北京航空航天大学2020届研究生《数值分析》实验作业第九题院系：xx学院学号：姓名：2020年11月Q9:方程组A.4一、算法设计方案（一）总体思路1.题目要求∑∑===k i kj s r rsy x cy x p 00),(对f(x, y) 进行拟合，可选用乘积型最小二乘拟合。

),(i i y x 与),(i i y x f 的数表由方程组与表A-1得到。

2.),(**j i y x f 与1使用相同方法求得，),(**j i y x p 由计算得出的p(x,y)直接带入),(**j i y x 求得。

1. ),(i i y x 与),(i i y x f 的数表的获得对区域D ={ （x,y)|1≤x ≤1.24,1.0≤y ≤1.16}上的f (x , y )值可通过xi=1+0.008i ，yj=1+0.008j ，得到),(i i y x 共31×21组。

将每组带入A4方程组，即可获得五个二元函数组，通过简单牛顿迭代法求解这五个二元数组可获得z1~z5有关x,y 的表达式。

再将),(i i y x 分别带入z1~z5表达式即可获得f(x,y)值。

2.乘积型最小二乘曲面拟合2.1使用乘积型最小二乘拟合，根据k 值不用，有基函数矩阵如下：⎪⎪⎪⎭⎫ ⎝⎛=k i i k x x x x B 0000 ， ⎪⎪⎪⎭⎫ ⎝⎛=k j jk y y y y G 0000数表矩阵如下：⎪⎪⎪⎭⎫⎝⎛=),(),(),(),(0000j i i j y x f y x f y x f y x f U记C=[rs c ]，则系数rs c 的表达式矩阵为：11-)(-=G G UG B B B C T TT ）（通过求解如下线性方程，即可得到系数矩阵C 。

UG B G G C B B T T T =)()(2.2计算),(),,(****j i j i y x p y x f (i =1,2,…,31 ; j =1,2,…,21) 的值),(**j i y x f 的计算与),(j i y x f 相同。

北航数值分析大作业3(学硕)

《数值分析》作业三院系：机械学院学号：SY1307145姓名：龙安林2013年11 月24 日1. 算法设计1) 开始；2) 计算数组[][]0.08,0.050.5,0,1,2,,10;0,1,2,,20x i i y j j i j ==+=⋯=⋯（）； 3) 将点[][],0,1,2,,10;0,1,2,,20x i y j i j =⋯=⋯（）,（）带入非线性方程组： 0.5cos 2.670.5sin 1.070.5cos 3.740.5sin 0.79t u v w x t u v w y t u v w x t u v w y +++-=⎧⎪+++-=⎪⎨+++-=⎪⎪+++-=⎩ 得出相应的点,t u （）； 4) 选择拉格朗日插值法，将,t u （）作为中间变量，在题目所给出的二维数表中进行二次代数插值，得到[][],)(z f x i y j =；5) 输出数表：[][][][]()()0,1,2,,10;0,1,2,,20,,,x i y j f x i y j i j =⋯=⋯； 6) 令k=0；7) 以()()(),,,0,1,r r r s x x y y r s ϕψ===…,k 为拟合基函数，将上述数表作为拟合条件，对于给定的k 值，得到矩阵B 、G 、U ；8) 令-1-1(),()T T T A B B B U C AG G G ==，用选主元的LU 分解法分别计算矩阵A 和C 的各列，最后得到系数矩阵C ；9) 以公式：()()()00,k ki j rs r i s j s r p x y C x y ϕψ===∑∑计算每个点的拟合值；10) 利用公式：()()()2102000,,i j i j i j f x y p x y σ===-∑∑计算拟合误差，当σ≤10-7时，循环结束，否则k=k+1，转（6）；11) 令[][]()**0.10.50.2 1,2,81,2,5x i i y j j i j ==+=⋯=⋯；，；，；12) 计算()()()******,,,,,i j i j i jf x y p x y delta x y ，输出数表，观察逼近效果； 13) 结束。

北航研究生数值分析大作业三

数值分析—计算实习作业三学院：17系专业：精密仪器及机械姓名：张大军学号：DY1417114一、程序设计方案程序设计方案流程图如图1所示。

（注：由本人独立完成，并且有几处算法很巧妙，但同时也有许多不足，可以优化和模块化，由于时间原因只实现了调试通过）图1.程序设计方案流程图二、程序源代码#include <iostream.h>#include <iomanip.h>#include <math.h>#include<stdio.h>#include <conio.h>#define M 10000#define N 4#define E 1.0e-12int zuixiaci;static double c[9][9];static double bijin[8][5];int main(){double X[N]={0,0,0,1};double T[11][21],U[11][21],xianshi[11][21];double diertX[N];double F[N];double f[N][N];double Max1=0,Max2=0;int k,i,j,t,tt=0,yao=0;void qiuF(double * X,double *F,int i,int j);void qiuF2(double *X,double *F,int i,int j);void qiuf(double * X,double (*f)[N]);void qiudiertX(double (*a)[N],double*b,double*X); double gouzaohs(double t,double u); void solve_C(double (*U)[21]); void pp(double (*U)[21],int k);for(i=0;i<11;i++)for(j=0;j<21;j++){for(k=0;k<M;k++){qiuF(X,F,i,j);qiuf(X,f);qiudiertX(f,F,diertX);for(t=0;t<N;t++){X[t]=X[t]+diertX[t];}Max1=0,Max2=0;for(t=0;t<N;t++){if(Max1<fabs(X[t]))Max1=fabs(X[t]);if(Max2<fabs(diertX[t]))Max2=fabs(diertX[t]);}if((Max2/Max1)<=E){k=M;yao=1;T[i][j]=X[0];U[i][j]=X[1];xianshi[i][j]=gouzaohs(X[0],X[1]);cout<<setiosflags(ios::scientific)<<setprecision(12);cout<<setprecision(2)<<"("<<setw(5)<<0.08*i<<","<<setw(5)<< (0.5+0.05*j)<<",";cout<<setprecision(12)<<setw(21)<<xianshi[i][j]<<") ";if(tt==3){tt=0;cout<<'\n';cout<<'\n';}else{tt++;}}}if(yao==0)cout<<"迭代不成功"<<endl; yao=0;}cout<<endl;solve_C(xianshi);pp(xianshi,zuixiaci);tt=0;for(i=1;i<9;i++)for(j=1;j<6;j++){for(k=0;k<M;k++){qiuF2(X,F,i,j);qiuf(X,f);qiudiertX(f,F,diertX);for(t=0;t<N;t++){X[t]=X[t]+diertX[t];}Max1=0,Max2=0;for(t=0;t<N;t++){if(Max1<fabs(X[t]))Max1=fabs(X[t]);if(Max2<fabs(diertX[t]))Max2=fabs(diertX[t]);}if((Max2/Max1)<=E){k=M;yao=1;xianshi[i-1][j-1]=gouzaohs(X[0],X[1]);cout<<setiosflags(ios::scientific)<<setprecision(12);cout<<setprecision(2)<<"("<<setw(5)<<0.1*i<<","<<setw(5)<<( 0.5+0.2*j)<<",";cout<<setprecision(12)<<setw(21)<<xianshi[i-1][j-1]<<","<<set w(21)<<bijin[i-1][j-1]<<") ";if(tt==2){tt=0;cout<<'\n';}else{tt++;}}}if(yao==0)cout<<"迭代不成功"<<endl;yao=0;}cout<<endl;return 1;}void qiuF(double *X,double *F,int i,int j){F[0]=-(0.5*cos(X[0])+X[1]+X[2]+X[3]-0.08*i-2.67);F[1]=-(X[0]+0.5*sin(X[1])+X[2]+X[3]-(0.5+0.05*j)-1.07);F[2]=-(0.5*X[0]+X[1]+cos(X[2])+X[3]-0.08*i-3.74);F[3]=-(X[0]+0.5*X[1]+X[2]+sin(X[3])-(0.5+0.05*j)-0.79); }void qiuF2(double *X,double *F,int i,int j){F[0]=-(0.5*cos(X[0])+X[1]+X[2]+X[3]-0.1*i-2.67);F[1]=-(X[0]+0.5*sin(X[1])+X[2]+X[3]-(0.5+0.2*j)-1.07);F[2]=-(0.5*X[0]+X[1]+cos(X[2])+X[3]-0.1*i-3.74);F[3]=-(X[0]+0.5*X[1]+X[2]+sin(X[3])-(0.5+0.2*j)-0.79); }void qiuf(double *X,double (*f)[N]){f[0][0]=-0.5*sin(X[0]);f[0][1]=1;f[0][2]=1;f[0][3]=1;f[1][0]=1;f[1][1]=0.5*cos(X[1]);f[1][2]=1;f[1][3]=1;f[2][0]=0.5;f[2][1]=1;f[2][2]=-sin(X[2]);f[2][3]=1;f[3][0]=1;f[3][1]=0.5;f[3][2]=1;f[3][3]=cos(X[3]);}//求解关于变化X的线性方程组void qiudiertX(double (*a)[N],double*b,double*X) {double H[N][N]={0},l[N]={0};double B;double sum;int i,j,m,k,z;for(k=0;k<N-1;k++){for(j=k;j<N;j++){l[j]=a[k][j];}z=k;for(m=k;m<N;m++){if(fabs(a[z][k])<fabs(a[m][k]))z=m;}for(j=k;j<N;j++){a[k][j]=a[z][j];a[z][j]=l[j];}B=b[k];b[k]=b[z];b[z]=B;for(i=k+1;i<N;i++){H[i][k]=a[i][k]/a[k][k];for(j=k+1;j<N;j++)a[i][j]=a[i][j]-H[i][k]*a[k][j];b[i]=b[i]-H[i][k]*b[k];}}if(a[N-1][N-1]==0){cout<<"算法失效，停止计算"<<endl; }else{X[N-1]=b[N-1]/a[N-1][N-1];for(k=N-2;k>=0;k--){sum=0;for(j=k+1;j<N;j++){sum=sum+a[k][j]*X[j];}X[k]=(b[k]-sum)/a[k][k];}}}//作二元差值，使用分片二次代数插值double gouzaohs(double t,double u){double T[6]={0,0.2,0.4,0.6,0.8,1},U[6]={0,0.4,0.8,1.2,1.6,2};double Z[6][6]={-0.5,-0.34,0.14,0.94,2.06,3.5,-0.42,-0.5,-0.26,0.3,1.18,2.38,-0.18,-0.5,-0.5,-0.18,0.46,1.42,0.22,-0.34,-0.58,-0.5,-0.1,0.62,0.78,-0.02,-0.5,-0.66,-0.5,-0.02,1.5,0.46,-0.26,-0.66,-0.74,-0.5};double g=0,sum=0,sum1=1,sum2=1;int i=0,j=0,k=0,r=0,kk=0,rr=0;for(i=1;(i<6)&&(T[i]-0.1<t);i++){}for(j=1;(j<6)&&(U[j]-0.2<u);j++){}if(i==1)i=2;if(i==6)i=5;if(j==1)j=2;if(j==6)j=5;sum=0;for(k=i-2;k<i+1;k++)for(r=j-2;r<j+1;r++){sum1=1;sum2=1;for(kk=i-2;kk<i+1;kk++){if(k!=kk){sum1=sum1*(t-T[kk])/(T[k]-T[kk]);}}for(rr=j-2;rr<j+1;rr++){if(r!=rr){sum2=sum2*(u-U[rr])/(U[r]-U[rr]);}}sum=sum+sum1*sum2*Z[k][r];}g=sum;return g;}//求r*s阶矩阵A与s*t阶矩阵B的乘积矩阵Cvoid Multi(double *a, double *b, double *c, int la, int lb, int lc, int r, int s, int t){int i, j, k;for (i=0; i<r; i++)for (j=0; j<t; j++){*(c+i*lc+j)=0;for (k=0; k<s; k++)*(c+i*lc+j)+=*(a+i*la+k)*(*(b+k*lb+j));}}//求n阶方阵A的逆矩阵Bdouble Inverse(double *a, double *b, int la, int lb, int n){int i, j, k;double temp;for(i=0; i<n; i++)for(j=0; j<n; j++)if (i==j)*(b+i*lb+j)=1;else*(b+i*lb+j)=0;for (k=0; k<n; k++){j=k;for (i=k+1; i<n; i++)if (fabs(*(a+i*la+k))>fabs(*(a+j*la+k))) j=i;if (j!=k)for (i=0; i<n; i++){temp=*(a+j*la+i);*(a+j*la+i)=*(a+k*la+i);*(a+k*la+i)=temp;temp=*(b+j*lb+i);*(b+j*lb+i)=*(b+k*lb+i);*(b+k*lb+i)=temp;}if (*(a+k*la+k)==0)return 0;if ((temp=*(a+k*la+k))!=1)for (i=0; i<n; i++){*(a+k*la+i)/=temp;*(b+k*lb+i)/=temp;}for (i=0; i<n; i++)if ((*(a+i*la+k)!=0) && (i!=k)){temp=*(a+i*la+k);for (j=0; j<n; j++){*(a+i*la+j)-=temp*(*(a+k*la+j));*(b+i*lb+j)-=temp*(*(b+k*lb+j));}}}return 0;}void solve_C(double (*U)[21]){int i,j,r,s,k;double t1[21][21], t2[21][21], t3[21][21],d[9][9],e[9][9];double B[11][9], B_T[9][11], G[21][9], G_T[9][21],P[11][21];double temp, FangCha;for(i=0;i<9;i++){for(j=0;j<11;j++){B[j][i]=pow(0.08*j,i);B_T[i][j]=pow(0.08*j,i);}for(j=0;j<21;j++){G[j][i]=pow(0.5+0.05*j,i);G_T[i][j]=pow(0.5+0.05*j,i);}}for (k=0; k<9; k++){FangCha=0;Multi(B_T[0], B[0], t1[0], 11, 9, 21, k+1, 11, k+1);Inverse(t1[0], c[0], 21, 9, k+1);Multi(e[0], c[0], d[0], 9, 9, 9, k+1, k+1, k+1);Multi(c[0], B_T[0], t1[0], 9, 11, 21, k+1, k+1, 11);Multi(t1[0], U[0], t2[0], 21, 21, 21, k+1, 11, 21);Multi(G_T[0], G[0], t1[0], 21, 9, 21, k+1, 21, k+1);Inverse(t1[0], c[0], 21, 9, k+1);Multi(G[0], c[0], t3[0], 9, 9, 21, 21, k+1, k+1);Multi(t2[0], t3[0], c[0], 21, 21, 9, k+1, 21, k+1);for(i=0;i<11;i++)for(j=0;j<21;j++){temp=0;for(r=0;r<k+1;r++)for(s=0;s<k+1;s++)temp+=c[r][s]*B[i][r]*G[j][s];P[i][j]=temp;FangCha+=(U[i][j]-temp)*(U[i][j]-temp);}cout<<"k="<<setw(5)<<k<<";"<<setw(5)<<"Sigma="<<FangCha<<" ;\n"<<'\n';if(FangCha<=1.0e-7){zuixiaci=k;cout<<"达到精度要求时: k="<<setw(5)<<k<<";"<<setw(5)<<"Sigma="<<FangCha<<";\n";cout<<" 系数c(r,s)如下：\n";for(i=0;i<k+1;i++){for(j=0;j<k+1;j++){cout<<"C("<<i<<","<<j<<")="<<setw(21)<<c[i][j]<<"; ";}cout<<endl<<'\n';}cout<<endl;return;}}cout<<"经过8次拟合没有达到所需精度；"<<endl;//最高可拟合10次return;}void pp(double (*U)[21],int k){int i,j,r,s;double B[8][9],G[5][9],temp;for(i=0;i<k+1;i++){for(j=0;j<8;j++){B[j][i]=pow(0.1*(j+1),i);}for(j=0;j<5;j++){G[j][i]=pow(0.5+0.2*(j+1),i);}}for(i=0;i<8;i++)for(j=0;j<5;j++){temp=0;for(r=0;r<k+1;r++)for(s=0;s<k+1;s++)temp+=c[r][s]*B[i][r]*G[j][s];bijin[i][j]=temp;}}三、程序运行结果显示程序运行结果显示如图2。

北航航空发动机原理3大作业

航空发动机原理Ⅲ大作业—发动机设计点热力计算学院能源与动力工程学院一. 设计要求1.完成一台发动机的设计点热力计算1）完成发动机循环参数的选取2）完成发动机各部件设计参数（包括冷却空气量及其分配关系）的选取3）说明以上参数选取的具体理由和依据4）完成发动机各部件进出口截面参数（流量总）完成发动机各部件进出口截面参数（流量、总温、总压）的计算5）完成发动机总性能（推力、耗油率）的计算，并满足给定的要求（误差并满足给定的要求（误差±2%）2.题目：分排涡扇发动机，高度11km，马赫数0.8，标准大气条件下，发动机推力2500daN，耗油率耗油率0.6kg/（daN.h）二.设计参数1. 设计点参数设计点物性参数空气比热Cp：1.005KJ/Kg燃气比热Cpg：1.244KJ/Kg空气绝热指数k：1.42.发动机参数（资料参考）3.设计点飞行条件4.部件效率和损失系数高压轴机械效率：ηmH=0.98低压轴机械效率：ηmL=0.98高压涡轮相对冷气量：δ1=7%低压涡轮相对冷气量：δ2=1%飞机引气量：β=1%相对功率提取效率：相对功率提取系数：CT0=3三.循环参数的初步选取范围1.涵道比随着涵道比B的增加，当单位推力一定时，存在最佳涵道比，使sfc达到最小值，而T t4随涵道比单调增加，因此B过大或者过小会使sfc达不到要求，且B过大会使涡轮前温度超温，当单位推力较小时，sfc随B的变化曲线在附近较为平坦，因此减小B，并不严重增加sfc，但可使涡轮前总温T t4显著降低。

根据资料查得的发动机参数，初始可取涵道比B=6~12。

2.涡轮前温度根据现有涡轮材料和冷却技术水平，涡轮前温度最高能达到2200K，且在亚声速飞行时，涡轮前温度过高会使耗油率增加。

根据现有发动机参数，选取涡轮前温度。

3.风扇增压比风扇增压比一般随涵道比增加而降低，对于涵道比为B=6~10的涡扇发动机，一般取。

4. 总增压比π在给定涡轮前温度前提下，存在使单位推力达到最大值的最佳增压比，且随涡轮前温度提高而增大；存在使耗油率达到最小值的压气机最经济增压比。

北航数值分析大作业第一题幂法与反幂法

北航数值分析大作业第一题幂法与反幂法-标准化文件发布号：（9456-EUATWK-MWUB-WUNN-INNUL-DDQTY-KII《数值分析》计算实习题目第一题：1. 算法设计方案（1）1λ，501λ和s λ的值。

1)首先通过幂法求出按模最大的特征值λt1，然后根据λt1进行原点平移求出另一特征值λt2，比较两值大小，数值小的为所求最小特征值λ1，数值大的为是所求最大特征值λ501。

2)使用反幂法求λs ，其中需要解线性方程组。

因为A 为带状线性方程组，此处采用LU 分解法解带状方程组。

（2）与140k λλμλ-5011=+k 最接近的特征值λik 。

通过带有原点平移的反幂法求出与数k μ最接近的特征值 λik 。

（3）2cond(A)和det A 。

1）1=nλλ2cond(A)，其中1λ和n λ分别是按模最大和最小特征值。

2）利用步骤（1）中分解矩阵A 得出的LU 矩阵，L 为单位下三角阵,U 为上三角阵，其中U 矩阵的主对角线元素之积即为det A 。

由于A 的元素零元素较多，为节省储存量，将A 的元素存为6×501的数组中，程序中采用get_an_element()函数来从小数组中取出A 中的元素。

2.全部源程序#include <stdio.h>#include <math.h>void init_a();//初始化Adouble get_an_element(int,int);//取A 中的元素函数double powermethod(double);//原点平移的幂法double inversepowermethod(double);//原点平移的反幂法int presolve(double);//三角LU分解int solve(double [],double []);//解方程组int max(int,int);int min(int,int);double (*u)[502]=new double[502][502];//上三角U数组double (*l)[502]=new double[502][502];//单位下三角L数组double a[6][502];//矩阵Aint main(){int i,k;double lambdat1,lambdat2,lambda1,lambda501,lambdas,mu[40],det;double lambda[40];init_a();//初始化Alambdat1=powermethod(0);lambdat2=powermethod(lambdat1);lambda1=lambdat1<lambdat2?lambdat1:lambdat2;lambda501=lambdat1>lambdat2?lambdat1:lambdat2;presolve(0);lambdas=inversepowermethod(0);det=1;for(i=1;i<=501;i++)det=det*u[i][i];for (k=1;k<=39;k++){mu[k]=lambda1+k*(lambda501-lambda1)/40;presolve(mu[k]);lambda[k]=inversepowermethod(mu[k]);}printf("------------所有特征值如下------------\n");printf("λ=%1.11e λ=%1.11e\n",lambda1,lambda501);printf("λs=%1.11e\n",lambdas);printf("cond(A)=%1.11e\n",fabs(lambdat1/lambdas));printf("detA=%1.11e \n",det);for (k=1;k<=39;k++){printf("λi%d=%1.11e ",k,lambda[k]);if(k % 3==0) printf("\n");}delete []u;delete []l;//释放堆内存return 0;}void init_a()//初始化A{int i;for (i=3;i<=501;i++) a[1][i]=a[5][502-i]=-0.064;for (i=2;i<=501;i++) a[2][i]=a[4][502-i]=0.16;for (i=1;i<=501;i++) a[3][i]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i); }double get_an_element(int i,int j)//从A中节省存储量的提取元素方法{if (fabs(i-j)<=2) return a[i-j+3][j];else return 0;}double powermethod(double offset)//幂法{int i,x1;double u[502],y[502];double beta=0,prebeta=-1000,yita=0;for (i=1;i<=501;i++)u[i]=1,y[i]=0;//设置初始向量u[]for (int k=1;k<=10000;k++){yita=0;for (i=1;i<=501;i++) yita=sqrt(yita*yita+u[i]*u[i]);for (i=1;i<=501;i++) y[i]=u[i]/yita;for (x1=1;x1<=501;x1++){u[x1]=0;for (int x2=1;x2<=501;x2++)u[x1]=u[x1]+((x1==x2)(get_an_element(x1,x2)-offset):get_an_element(x1,x2))*y[x2];}prebeta=beta;beta=0;for (i=1;i<=501;i++) beta=beta+ y[i]*u[i];if (fabs((prebeta-beta)/beta)<=1e-12) {printf("offset=%f lambda=%f err=%e k=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};//输出中间过程,包括偏移量,误差,迭代次数}return (beta+offset);}double inversepowermethod(double offset)//反幂法{int i;double u[502],y[502];double beta=0,prebeta=0,yita=0;for (i=1;i<=501;i++)u[i]=1,y[i]=0; //设置初始向量u[]for (int k=1;k<=10000;k++){yita=0;for (i=1;i<=501;i++) yita=sqrt(yita*yita+u[i]*u[i]);for (i=1;i<=501;i++) y[i]=u[i]/yita;solve(u,y);prebeta=beta;beta=0;for (i=1;i<=501;i++) beta=beta+ y[i]*u[i];beta=1/beta;if (fabs((prebeta-beta)/beta)<=1e-12) {printf("offset=%f lambda=%f err=%e k=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};//输出中间过程,包括偏移量,误差,迭代次数}return (beta+offset);}int presolve(double offset)//三角LU分解{int i,k,j,t;double sum;for (k=1;k<=501;k++)for (j=1;j<=501;j++){u[k][j]=l[k][j]=0;if (k==j) l[k][j]=1;} //初始化LU矩阵for (k=1;k<=501;k++){for (j=k;j<=min(k+2,501);j++){sum=0;for (t=max(1,max(k-2,j-2)) ; t<=(k-1) ; t++)sum=sum+l[k][t]*u[t][j];u[k][j]=((k==j)(get_an_element(k,j)-offset):get_an_element(k,j))-sum;}if (k==501) continue;for (i=k+1;i<=min(k+2,501);i++){sum=0;for (t=max(1,max(i-2,k-2));t<=(k-1);t++)sum=sum+l[i][t]*u[t][k];l[i][k]=(((i==k)(get_an_element(i,k)-offset):get_an_element(i,k))-sum)/u[k][k];}}return 0;}int solve(double x[],double b[])//解方程组{int i,t;double y[502];double sum;y[1]=b[1];for (i=2;i<=501;i++){sum=0;for (t=max(1,i-2);t<=i-1;t++)sum=sum+l[i][t]*y[t];y[i]=b[i]-sum;}x[501]=y[501]/u[501][501];for (i=500;i>=1;i--){sum=0;for (t=i+1;t<=min(i+2,501);t++)sum=sum+u[i][t]*x[t];x[i]=(y[i]-sum)/u[i][i];}return 0;}int max(int x,int y){return (x>y?x:y);}int min(int x,int y){return (x<y?x:y);}3.计算结果结果如下图所示:部分中间结果：给出了偏移量(offset)，误差(err)，迭代次数(k)4.讨论迭代初始向量的选取对计算结果的影响,并说明原因使用u[i]=1（i=1,2,...,501）作为初始向量进行迭代，可得出以上结果。

北航数值分析实习题目第一题

《数值分析B》大作业一ZY1515105 樊雪松一．算法设计方案：1.矩阵A的存储与检索将带状线性矩阵A[501][501]转存为一个矩阵MatrixC[5][501] 。

在数组MatrixC[5][501]中检索A的带内元素a ij的方法是：A的带内元素a ij=C中的元素c i-j+2,j。

2.求解λ1，λ501，λs1、首先分别使用幂法和反幂法迭代求出矩阵按摸最大和最小的特征值λmax和λmin。

λmin即为λs；如果λ max>0,则λ501=λmax；如果λmax<0,则λ1=λmax。

2、使用带原点平移的幂法（mifa（）函数），令平移量p=λmax，求出对应的按摸最大的特征值λ’max，如果λ max>0,则λ1=λ’max+p；如果λmax<0,则λ501=λ’max+p。

3、求解A的与数μk=λ1+k（λ501-λ1）/40 的最接近的特征值λik （k=1,2，…，39）。

使用带原点平移的反幂法，令平移量p=μk，即可求出与μk最接近的特征值λ ik。

4、求解A的（谱范数）条件数cond（A）2和行列式detA。

cond（A）2=|λ1/λn|，其中λ1和λn分别是矩阵A的模最大和最小特征值。

求解矩阵A的行列式，可先对矩阵A进行LU分解后，detA等于U所有对角线上元素的乘积。

二．源程序#include<stdio.h>#include<math.h>#include<conio.h>//定义A中元素double C[5][501];double a[501];double b;double c;//声明所有函数void YaSuoJZ(double C[5][501],double a[501],double b,double c) ;//压缩矩阵函数double mifa(double C[5][501]); //幂法函数void daizhuangLU(double A[5][501]); //带状矩阵的LU分解double fanmifa(double C[5][501]);//反幂法函数//最值函数int max2(int x,int y);int max3(int x,int y,int z);int min(int x,int y);//最值函数int max2(int x,int y) //求2个数的最大值{int z;z=x>y?x:y;return(z);}int max3(int x,int y,int z) //求3个数的最大值{int w;w = z > max2(x,y)? z:max2(x,y);return(w);}int min(int x,int y) //求2个数的最小值{int z;z=x>y?y:x;return(z);}//将矩阵A压缩存储在矩阵C中void YaSuoJZ(double C[5][501],double a[501],double b,double c) {int i;for(i=0;i<=500;i++){if(i>=2) C[0][i]=c;else C[0][i]=0;if(i>=1) C[1][i]=b;else C[1][i]=0;if(i<=499) C[3][i]=b;else C[3][i]=0;if(i<=498) C[4][i]=c;else C[4][i]=0;C[2][i]=a[i];}}//幂法函数：用幂法求矩阵模最大的特征值double mifa(double C[5][501]){double u[501];double y[501]={0},η=0;double β,βk=0;double ε=1;// ε为精度double sumu=0,sumAY=0;int i,j,k=1; //k为循环次数for (i=0;i<=500;i++) //取任一非零向量u0u[i] = 1.0;while(ε>=1e-12){for(i=0;i<=500;i++) //求u(k-1)的2范数ηsumu=sumu+u[i]*u[i];η=sqrt(sumu);sumu=0;for(i=0;i<=500;i++) //求y(k-1)y[i]=u[i]/η;for(i=0;i<=500;i++) //求u(k)的各分量u[i]{for(j=max2(0,i-2);j<=min(i+2,500);j++)sumAY=sumAY+C[i-j+2][j]*y[j];u[i]=sumAY;sumAY=0;}//求幂法中的βkβ=βk; //将β(k-1)放在β中βk=0;for(i=0;i<=500;i++) //求βkβk=βk+y[i]*u[i];if(k>=2)ε=fabs(βk-β)/fabs(βk);k++;}return(βk);}//带状矩阵的LU分解void daizhuangLU(double A[5][501]){int i,j,k,m,t;double sumukj=0,sumlik=0;for(k=0;k<=500;k++){for(j=k;j<=min(k+2,500);j++) //求ukj并存在A[k-j+2][j]中{for(t=max3(0,k-2,j-2);t<=k-1;t++)sumukj=sumukj+A[k-t+2][t]*A[t-j+2][j];A[k-j+2][j]=A[k-j+2][j]-sumukj;sumukj=0;}if(k<500)for(i=k+1;i<=min(k+2,500);i++) //求lik并存在A[i-k+2][k]中{for(m=max3(0,i-2,k-2);m<=k-1;m++)sumlik=sumlik+A[i-m+2][m]*A[m-k+2][k];A[i-k+2][k]=(A[i-k+2][k]-sumlik)/A[2][k];sumlik=0;}}}//反幂法函数：用反幂法求矩阵的模最小的特征值double fanmifa(double M[5][501]){double u[501];double y[501]={0},x[501],η=0;double fβ,fβk=0;double ε=1;double fsumu=0,sumLX=0,sumUu=0;int i,t,m,k=1;for(i=0;i<=500;i++) //任取一非零向量u0u[i]=1;daizhuangLU(M); //对A进行LU分解A=LU，Au(k)=y(k-1)等价于Uu(k)=x和Lx=y(k-1) while(ε>=1e-12){for(i=0;i<=500;i++) //求u(k-1)的2范数ηfsumu=fsumu+u[i]*u[i];η=sqrt(fsumu);fsumu=0;for(i=0;i<=500;i++) //求y(k-1)y[i]=u[i]/η;for(i=0;i<=500;i++) //求中间向量xx[i]=y[i];for(i=1;i<=500;i++){for(t=max2(0,i-2);t<=i-1;t++)sumLX=sumLX+M[i-t+2][t]*x[t];x[i]=x[i]-sumLX;sumLX=0;}u[500]=x[500]/C[2][500]; //求u(k)的各分量u[i]for(i=499;i>=0;i--){for(m=i+1;m<=min(i+2,500);m++)sumUu=sumUu+M[i-m+2][m]*u[m];u[i]=(x[i]-sumUu)/M[2][i];sumUu=0;}//求反幂法中的βkfβ=fβk; //将fβ(k-1)放在fβ中fβk=0;for(i=0;i<=500;i++) //求fβkfβk=fβk+y[i]*u[i];if(k>=2)ε=fabs(1/fβk-1/fβ)/fabs(1/fβk);k++;}return(1/fβk);}//主函数void main(){int i,j,k;double λ1,λ501,λm,λm1,λm2,λs,λ,p;double cond,detA=1;for(i=1;i<=501;i++)a[i-1]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i);b=0.16;c=-0.064;YaSuoJZ(C,a,b, c); //将矩阵A中元素压缩存储在C中λm1=mifa(C); //对A用幂法求出模最大的特征值λm1λs=fanmifa(C); //对A用反幂法求出模最小的特征值λsYaSuoJZ(C,a,b, c); //还原矩阵A中元素并压缩存储在C中for(j=0;j<=500;j++) //对A进行平移，平移量为λm1，平移后矩阵元素压缩存储在C中C[2][j]=C[2][j]-λ?m1;λm=mifa(C);λm2=λm1+λm; //λm1与λm2是矩阵的最大最小特征值if(λm1>λm2) //判断A最大最小特征值{λ501=λm1;λ1=λm2;}else{λ501=λm2;λ1=λm1;}printf("数值分析计算实习第一题\n\n ZY1515105 樊雪松\n\n （1）A的最大最小以及模最小的特征值\n");printf("A的最小特征值λ1=%.13e\n",λ1);printf("A的最大特征值λ501=%.13e\n",λ501);printf("A的模最小特征值λs=%.13e\n",λs);printf("\n（2）与数μk最接近的特征值\n");printf("\t要求接近的值\t\t\t实际求得的特征值\n");YaSuoJZ(C,a,b, c); //还原矩阵A中元素并压缩存储在C中for(k=1;k<=39;k++){p=λ1+k*(λ501-λ1)/40;for(j=0;j<=501;j++)C[2][j]=C[2][j]-p;λ=fanmifa(C)+p;printf("μ%d=%.13e λ%d=%.13e\n",k,p,k,λ);YaSuoJZ(C,a,b, c); //还原矩阵A中元素并压缩存储在C中}printf("\n（3）计算A的条件数cond(A)和行列式detA\n");cond=λm1/λs;daizhuangLU(C);for(j=0;j<=500;j++)detA=detA*C[2][j];printf("A的条件数cond(A)=%.13e\n",cond);printf("A的行列式detA=%.13e\n",detA);getch();}三、运行结果数值分析计算实习第一题ZY1515105 樊雪松（1）A的最大最小以及模最小的特征值A的最小特征值λ1=-1.0700113615018e+001A的最大特征值λ501=9.7246340987773e+000A的模最小特征值λs=-5.5579107942295e-003（2）与数μk最接近的特征值要求接近的值实际求得的特征值μ1=-1.0189494922173e+001 λ1=-1.0182934033146e+001 μ2=-9.6788762293280e+000 λ2=-9.5857074250676e+000 μ3=-9.1682575364831e+000 λ3=-9.1726724239280e+000 μ4=-8.6576388436383e+000 λ4=-8.6522840078976e+000 μ5=-8.1470201507934e+000 λ5=-8.0934838086753e+000 μ6=-7.6364014579485e+000 λ6=-7.6594054076924e+000 μ7=-7.1257827651036e+000 λ7=-7.1196846486912e+000 μ8=-6.6151640722588e+000 λ8=-6.6117643393973e+000 μ9=-6.1045453794139e+000 λ9=-6.0661032265951e+000 μ10=-5.5939266865690e+000 λ10=-5.5851010526284e+000 μ11=-5.0833079937241e+000 λ11=-5.1140835298122e+000 μ12=-4.5726893008792e+000 λ12=-4.5788721768651e+000 μ13=-4.0620706080344e+000 λ13=-4.0964709262599e+000 μ14=-3.5514519151895e+000 λ14=-3.5542112157508e+000 μ15=-3.0408332223446e+000 λ15=-3.0410900181333e+000 μ16=-2.5302145294997e+000 λ16=-2.5339703111304e+000 μ17=-2.0195958366549e+000 λ17=-2.0032307695635e+000μ18=-1.5089771438100e+000 λ18=-1.5035576112274e+000μ19=-9.9835845096511e-001 λ19=-9.9355860600754e-001μ20=-4.8773975812023e-001 λ20=-4.8704267388496e-001μ21=2.2878934724645e-002 λ21=2.2317362495748e -002μ22=5.3349762756952e-001 λ22=5.3241747420686e -001μ23=1.0441163204144e+000 λ23=1.0528989626935e+000μ24=1.5547350132593e+000 λ24=1.5894458818809e+000μ25=2.0653537061042e+000 λ25=2.0603304602743e+000μ26=2.5759723989490e+000 λ26=2.5580755970728e+000μ27=3.0865910917939e+000 λ27=3.0802405093071e+000μ28=3.5972097846388e+000 λ28=3.6136208676923e+000μ29=4.1078284774837e+000 λ29=4.0913785104506e+000μ30=4.6184471703285e+000 λ30=4.6030353782791e+000μ31=5.1290658631734e+000 λ31=5.1329242838984e+000μ32=5.6396845560183e+000 λ32=5.5949063480833e+000μ33=6.1503032488632e+000 λ33=6.0809338570269e+000μ34=6.6609219417080e+000 λ34=6.6803540921116e+000μ35=7.1715406345529e+000 λ35=7.2938774481266e+000μ36=7.6821593273978e+000 λ36=7.7171117142356e+000μ37=8.1927780202427e+000 λ37=8.2252200140502e+000μ38=8.7033967130876e+000 λ38=8.6486660651935e+000μ39=9.2140154059324e+000 λ39=9.2542003445750e+000（3）计算A 的条件数cond(A)和行列式detAA 的条件数cond(A)=1.9252042739022e+003A 的行列式detA=2.7727861417521e+118四、结果分析设A 的n 个线性无关的特征向量为1x ，2x ,…,n x ，其相对应的特征值满足的关系为n λλλλ≥≥≥> 321。

北航数值分析大作业题目三

《数值分析》第三次大作业一、算法的设计方案：（一）、总体方案设计：（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 。

北航数值分析大作业第二题

“数值分析“计算实习大作业第二题——SY1415215孔维鹏一、计算说明本程序采用带双步位移的QR方法求解矩阵A的所有特征值，然后采用反幂法求解矩阵A的实特征值对应的特征向量。

在采用带双步位移的QR方法求解特征值时，对教材上所提供的具体算法作稍微的改动，以简化程序，具体算法如下所示：1、计算出A拟上三角化后的矩阵，给定精度水平和最大迭代次数L；2、记，令k=1，m=n；3、如果，则可直接计算出最后1或2个特征值，转8，否则转4；4、如果，则可得一个特征值，置m=m-1；转3，否则转5；5、如果，则可得两个特征值，置m=m-2；转3，否则转6；6、记，计算7、k=k+1,转38、A的全部特征值已经求出，停止计算。

二、计算源程序（FORTRAN）PROGRAM SY1415215_2PARAMETER (N=10)DIMENSION A(N,N),A1(N,N),A2(N,N),C(2,N),Q(N,N),R(N,N),CR(N),CM(N)!C为存储特征值的数组,1为实部，为虚部REAL(8) A,A1,A2,C,Q,R,CME=1E-12 !精度水平L=1000 !迭代最大次数OPEN(1,FILE='数值分析大作业第二题计算结果.TXT')DO I=1,NDO J=1,NIF(I==J) THENA(I,J)=*COS(I+*J)ELSEA(I,J)=SIN*I+*J)ENDIFENDDOENDDOA1=AWRITE(*,"('矩阵A为：')")WRITE(1,"('矩阵A为：')")DO I=1,NDO J=1,NWRITE(*,"(2X,,2X,\)") A(I,J)WRITE(1,"(2X,,2X,\)") A(I,J)ENDDOWRITE(*,"(' ')")WRITE(1,"(' ')")ENDDO!使用矩阵的拟上三角化的算法将矩阵A化为拟上三角矩阵A（n-1）CALL HESSENBERG(A,N)WRITE(*,"('拟上三角化后矩阵A（n-1）为：')")WRITE(1,"('拟上三角化后矩阵A（n-1）为：')")DO I=1,NDO J=1,NWRITE(*,"(2X,,2X,\)") A(I,J)WRITE(1,"(2X,,2X,\)") A(I,J)ENDDOWRITE(*,"('')")WRITE(1,"('')")ENDDO!计算对矩阵A（n-1）实行QR方法迭代结束后所得矩阵A2=ACALL QRD(A2,N,Q,R)WRITE(*,"('对矩阵A（n-1）实行QR方法迭代结束后所得Q为：')") WRITE(1,"('对矩阵A（n-1）实行QR方法迭代结束后所得Q为：')") DO I=1,NDO J=1,NWRITE(*,"(2X,,2X,\)") Q(I,J)WRITE(1,"(2X,,2X,\)") Q(I,J)ENDDOWRITE(*,"('')")WRITE(1,"('')")ENDDOWRITE(*,"('对矩阵A（n-1）实行QR方法迭代结束后所得R为：')") WRITE(1,"('对矩阵A（n-1）实行QR方法迭代结束后所得R为：')") DO I=1,NDO J=1,NWRITE(*,"(2X,,2X,\)") R(I,J)WRITE(1,"(2X,,2X,\)") R(I,J)ENDDOWRITE(*,"('')")WRITE(1,"('')")ENDDO!使用带双步位移的QR方法求解矩阵A（n-1）的特征值K=1M=NDO WHILE(K<=L)IF(M<=2) THENIF(M==1) THENC(1,M)=A(M,M)ELSE IF(M==2) THENCALL CALCUS(A,N,M,C)ENDIFEXITELSE IF(ABS(A(M,M-1))<E) THENC(1,M)=A(M,M)M=M-1ELSE IF(ABS(A(M-1,M-2))<E) THENCALL CALCUS(A,N,M,C)M=M-2ELSECALL CALM(A,M,N)ENDIFK=K+1ENDDOWRITE(*,"('矩阵A的全部特征值为：')")WRITE(1,"('矩阵A的全部特征值为：')")DO J=1,NWRITE(*,",'+',,'i')") C(1,J),C(2,J)WRITE(1,",'+',,'i')") C(1,J),C(2,J)ENDDO!使用反幂法求解A的相应于实特征值的特征向量J=1DO I=1,NIF(C(2,I)==0)THENCR(J)=C(1,I)J=J+1ENDIFENDDOJC=J-1WRITE(*,"('矩阵A的实特征值为：')")WRITE(1,"('矩阵A的实特征值为：')")DO I=1,JCWRITE(*,"") CR(I)WRITE(1,"") CR(I)ENDDODO II=1,JCDO I=1,NDO J=1,NIF(I==J) THENA(I,J)=A1(I,J)-CR(II)ELSEA(I,J)=A1(I,J)ENDIFENDDOENDDOCALL INPOVERMETHOD(A,N,CM)WRITE(*,"('与实特征值',,'对应的特征向量为：')") CR(II) WRITE(1,"('与实特征值',,'对应的特征向量为：')") CR(II) DO I=1,NWRITE(*,"(2X,,2X,\)") CM(I)WRITE(1,"(2X,,2X,\)") CM(I)ENDDOWRITE(*,"('')")WRITE(1,"('')")ENDDOCLOSE(1)END!***************拟上三角化子函数*************************！SUBROUTINE HESSENBERG(A,N)DIMENSION A(N,N),P(N),Q(N),W(N),U(N),AT(N,N)REAL(8) A,P,Q,W,U,ATREAL(8) S0,S1,S2,S3,S4,TDO L=1,N-2JUDGE=0DO I=L+2,NIF(A(I,L)/=0) THENJUDGE=1EXITENDIFENDDOIF(JUDGE==0) THENA=ACYCLEELSE IF(JUDGE/=0) THEN!计算DRS0=0DO I=L+1,NS0=S0+A(I,L)**2ENDDODR=SQRT(S0)!计算CRIF(A(L+1,L)==0)THENCR=DRELSECR=-SGN(A(L+1,L))*DRENDIF!计算HRHR=CR**2-CR*A(L+1,L)!给u赋值IF(I<L+1) THENU(I)=0ELSE IF(I==L+1) THEN U(I)=A(I,L)-CRELSE IF(I>L+1) THEN U(I)=A(I,L)ENDIFENDDO!计算PDO I=1,NDO J=1,NAT(I,J)=A(J,I)ENDDOENDDODO I=1,NS1=0DO J=1,NS1=S1+AT(I,J)*U(J)ENDDOP(I)=S1/HRENDDO!计算QDO I=1,NS2=0DO J=1,NS2=S2+A(I,J)*U(J)ENDDOQ(I)=S2/HRENDDO!计算TS3=0DO I=1,NS3=S3+P(I)*U(I)ENDDOT=S3/HR!计算WDO I=1,NW(I)=Q(I)-T*U(I)!计算A（r+1）DO I=1,NDO J=1,NA(I,J)=A(I,J)-W(I)*U(J)-U(I)*P(J)ENDDOENDDOENDIFENDDORETURNEND!***************符号函数子程序*****************！FUNCTION SGN(X)REAL(8) XIF(X>0) THENSGN=1ELSE IF(X<0) THENSGN=-1ELSE IF(X==0) THENSGN=0ENDIFEND!*********计算二阶子阵特征值s1，s2子函数*******！SUBROUTINE CALCUS(X,N,M,Y)DIMENSION X(N,N),Y(2,N)REAL(8) A,B,C,D,X,YA=1C=X(M-1,M-1)*X(M,M)-X(M-1,M)*X(M,M-1)B=-(X(M-1,M-1)+X(M,M))D=B**2-4*CIF(D>=0) THENY(1,M)=(-B-SQRT(D))/2Y(1,M-1)=(-B+SQRT(D))/2ELSEIF(D<0) THENY(1,M)=-B/2Y(1,M-1)=-B/2Y(2,M)=-SQRT(-D)/2Y(2,M-1)=-SQRT(-D)/2ENDIFRETURNEND!*********计算Mk，Ak+1子函数************！SUBROUTINE CALM(A,M,N)DIMENSION A(N,N),MK(M,M),X(M,M),QK(M,M),RK(M,M),S1(M,M),S2(M,M),QKT(M,M) REAL(8) A,MK,X,QK,RK,QKTREAL(8) S0,S1,S2DO I=1,MDO J=1,MIF(I==J) THENX(I,J)=1ELSEX(I,J)=0ENDIFENDDOENDDOS=A(M-1,M-1)+A(M,M)T=A(M-1,M-1)*A(M,M)-A(M,M-1)*A(M-1,M)DO I=1,MDO J=1,MS0=0DO K=1,MS0=S0+A(I,K)*A(K,J)ENDDOMK(I,J)=S0-S*A(I,J)+T*X(I,J)ENDDOENDDO!对Mk做QR分解CALL QRD(MK,M,QK,RK)DO I=1,MDO J=1,MQKT(I,J)=QK(J,I)ENDDOENDDODO I=1,MDO J=1,MS1(I,J)=0DO K=1,MS1(I,J)=S1(I,J)+QKT(I,K)*A(K,J)ENDDOENDDOENDDOA=S1DO I=1,MDO J=1,MS2(I,J)=0DO K=1,MS2(I,J)=S2(I,J)+A(I,K)*QK(K,J)ENDDOENDDOENDDOA=S2RETURNEND!************QR分解子程序***************！SUBROUTINE QRD(A,N,Q,R)DIMENSION A(N,N),AT(N,N),Q(N,N),U(N),W(N),P(N),R(N,N) REAL(8) A,AT,Q,U,W,P,RREAL(8) DR,S0,CR,HR,S1,S2DO I=1,NDO J=1,NIF(I==J) THENQ(I,J)=1ELSEQ(I,J)=0ENDIFENDDOENDDODO L=1,N-1JUDGE=0DO I=L+1,NIF(A(I,L)/=0) THENJUDGE=1EXITENDIF!A(I,L)中有一个不为零,判断条件为真，跳出循环转ENDDOIF(JUDGE==0) THENQ=QA=ACYCLE!A(I,L)全为零，结束本循环，进入下一个ELSE IF(JUDGE/=0) THEN!计算DRS0=0DO I=L,NS0=S0+A(I,L)**2ENDDODR=SQRT(S0)!计算CRIF(A(L,L)==0)THENCR=DRELSECR=-SGN(A(L,L))*DRENDIF!计算HRHR=CR**2-CR*A(L,L)!给u赋值DO I=1,NIF(I<L) THENU(I)=0ELSE IF(I==L) THENU(I)=A(I,L)-CRELSE IF(I>L) THENU(I)=A(I,L)ENDIFENDDO!计算WDO I=1,NS1=0DO J=1,NS1=S1+Q(I,J)*U(J)ENDDOW(I)=S1ENDDO!计算Q（r+1）DO I=1,NDO J=1,NQ(I,J)=Q(I,J)-W(I)*U(J)/HRENDDOENDDO!计算PDO I=1,NDO J=1,NAT(I,J)=A(J,I)ENDDOENDDODO I=1,NS2=0DO J=1,NS2=S2+AT(I,J)*U(J)ENDDOP(I)=S2/HRENDDO!计算A（r+1）DO I=1,NDO J=1,NA(I,J)=A(I,J)-U(I)*P(J)ENDDOENDDOENDIFENDDOQ=QR=ARETURNEND!*************运用反幂法求解矩阵A实特征值的特征向量***********！SUBROUTINE INPOVERMETHOD(A,N,Y)DIMENSION A(N,N),U(N),Y(N),U1(N,N),L1(N,N)REAL(8) E,Z,Z1,Z2,S1,S2,BREAL(8) A,U,Y,U1,L1U(I)=1ENDDO!任取非零向量UCALL DETA(A,N,U1,L1)Z2=EIZ1=E=1K=1DO WHILE (E>1E-12)S1=0DO I=1,NS1=S1+U(I)**2ENDDOB=SQRT(S1) !1DO I=1,NY(I)=U(I)/BENDDO!2CALL DOOLITTLE(U1,L1,Y,N,U) !3利用DOOLITTLE分解法法求解Au=yS2=0DO I=1,NS2=S2+Y(I)*U(I)ENDDOZ1=Z2Z2=S2 !4E=ABS(1/Z2-1/Z1)/ABS(1/Z2) !判断是否满足精度K=K+1ENDDORETURNENDSUBROUTINE DOOLITTLE(U,L,B1,N,X)DIMENSION B(N),U(N,N),X(N),Y(N),B1(N)REAL(8) L(N,N)REAL(8) B,U,X,Y,B1REAL(8) S1,S2,S3,S4B=B1Y(1)=B(1)S3=0DO M=1,I-1S3=S3+L(I,M)*Y(M)ENDDOY(I)=B(I)-S3ENDDOX(N)=Y(N)/U(N,N)DO I=N-1,1,-1S4=0DO M=I+1,NS4=S4+U(I,M)*X(M)ENDDOX(I)=(Y(I)-S4)/U(I,I)ENDDORETURNENDSUBROUTINE DETA(A1,N,U,L) DIMENSION A(N,N),U(N,N),A1(N,N) REAL(8) L(N,N)REAL(8) X,S1,S2REAL(8) A,U,A1X=1A=A1!对矩阵A进行Doolittle分解DO K=1,NDO J=K,NS1=0DO M=1,K-1S1=S1+L(K,M)*U(M,J)ENDDOU(K,J)=A(K,J)-S1A(K,J)=U(K,J)ENDDOIF (K==N) THENEXITELSEDO I=K+1,NS2=0DO M=1,K-1S2=S2+L(I,M)*U(M,K)ENDDOL(I,K)=(A(I,K)-S2)/U(K,K)A(I,K)=L(I,K)ENDDOENDIFENDDORETURNEND三、计算结果矩阵A为：+00 +00 +00 +00 +00 +00 +00 +00 +00 +00+00 +00 +00 +00 +00 +00 +00 +00 +00 +00+00 +00 +01 +00 +00 +00 +00 +00 +00+00 +00 +00 +01 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00+00 +00 +00 +00 +01 +00 +00 +00 +00+00 +00 +00 +00 +00 +00 +01 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00+00 +00 +00 +00 +00 +00 +00 +00 +01拟上三角化后矩阵A（n-1）为：+00 +01 +00 +00 +01 +00 +00 +00+01 +01 +01 +00 +00 +01 +00 +01 +00 +00+01 +01 +01 +00 +01 +00 +00 +00 +00+01 +01 +01 +01 +00 +00 +00 +00+01 +00 +00 +00 +00 +00+00 +01 +00 +00 +00 +00+00 +00 +00+00 +00 +00 +00+00 +00 +00+00 +00对矩阵A（n-1）实行QR方法迭代结束后所得Q为：+00 +00 +00 +00 +00 +00+00 +00 +00 +00+00 +00 +00 +00+00 +00 +00 +00 +00 +00+00 +00+00 +00 +00 +00 +00+00 +00 +00 +00+00+00 +00 +00+00 +00对矩阵A（n-1）实行QR方法迭代结束后所得R为：+01 +01 +01 +00 +01 +00 +00 +00 +00 +01 +00 +00 +00 +00 +00 +00 +00 +00 +01 +01 +00 +01 +00 +01 +00 +00 +00 +00 +01 +00 +00 +00 +00+00 +00 +01 +01 +00 +00 +00+00 +00 +01 +00 +00 +00 +00+00 +00 +00 +00 +00 +00 +00+00 +00 +00 +00 +01 +00 +00+00 +00 +00 +00 +00 +00+00 +00 +00 +00矩阵A的全部特征值为：+01+ +00i+01++00i+01++00i+01+ +00i+01+ +00i+00++00i+00++00i+00+ +00i+00+ +00i+ +00i矩阵A的实特征值为：+01+01+01+00+00与实特征值 +01对应的特征向量为：+00 +00 +00 +00 +00 +00 +00 +00 +00与实特征值 +01对应的特征向量为：+00 +00 +00 +00 +00 +00 +00与实特征值+01对应的特征向量为：+00 +00 +00与实特征值 +00对应的特征向量为：+00 +00 +00 +00 +00与实特征值 +00对应的特征向量为：+00 +00 +00 +00 +00 +00 +00与实特征值对应的特征向量为：+00 +00 +00 +00 +00 +00 +00 +00。

北航数值分析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 的效果。

北航数值分析全部三次大作业

北航数值分析大作业 第一题 幂法与反幂法

北航数值分析大作业3

北航研究生数值分析上机作业 三 (报告+所有程序大全)

数值分析大作业三四五六七

数值分析第三次作业

北航数值分析大作业三

北航数值分析报告大作业第三题(fortran)

BUAA数值分析大作业三

北航数值分析大作业3(学硕)

北航研究生数值分析大作业三

北航航空发动机原理3大作业

北航数值分析大作业第一题幂法与反幂法

北航数值分析实习题目第一题

北航数值分析大作业题目三

北航数值分析大作业第二题

北航数值分析A大作业3

北航数值分析大作业第一题幂法与反幂法

北航研究生数值分析上机作业三 (报告+所有程序大全)