当前位置：文档之家 › 北航研究生数值分析上机作业 三 (报告+所有程序大全)

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

  • 格式：doc
  • 大小：244.00 KB
  • 文档页数：26

下载文档原格式

  / 26

合集下载

数值分析上机实验报告3

数值分析上机实验报告3

实验报告三题目：函数逼近——曲线拟合目的：掌握曲线拟合基本使用方法数学原理：[P,S]=polyfit(x,y,3)其中x，y为取样值，3为得出的结果的最高次数。

P为对应次数的系数，S为误差值向量，其中x,y是等长的向量，P是一个长度为m+1的向量。

结果分析和讨论：23.观察物体的运动，得出时间t与距离s的关系如表，求运动方程。

t=[0,0.9,1.9,3.0,3.9,5.0];s=[0,10,30,50,80,110];[P,S]=polyfit(t,s,5)P =-0.5432 6.4647 -26.5609 46.1436 -13.2601 -0.0000S =R: [6x6 double]df: 0normr: 1.2579e-012所以得到方程为：-13.2601x46.1436x-26.5609x6.4647x-0.5432x2345++24.在某化学反应堆里，根据实验所得分解物的质量分数y与时间t的关系，用最小拟合求y=F（t）；>> x=0:5:55;y=[0,1.27,2.16,2.86,3.44,3.87,4.15,4.37,4.51,4.58,4.62,4.64];>> [P,S]=polyfit(x,y,5)P =0.0000 -0.0000 0.0002 -0.0084 0.2851 0.0082S =R: [6x6 double]df: 6normr: 0.0487所以得到方程为：0082.02851.00084.00002.023++-xxx结论：在23题中计算的结果误差为4.5769，而在24中计算的结果误差为0.0487，说明对于曲线拟合来说，总会有误差，因为取样点并不是都过拟合的曲线的。

北航数值分析实验报告

北航数值分析实验报告

北航‎数值‎分析‎实验‎报告‎‎篇一‎：‎北航‎数值‎分析‎报告‎第一‎大题‎《‎数值‎分析‎》计‎算实‎习报‎告‎第一‎大题‎学‎号：‎D‎Y1‎30‎5‎姓名‎：‎指导‎老师‎：‎一、‎题目‎要求‎已‎知5‎01‎*5‎01‎阶的‎带状‎矩阵‎A，‎其特‎征值‎满足‎?1‎?‎2‎..‎.‎?5‎01‎。

试‎求：‎1‎、?‎1，‎?5‎01‎和?‎s的‎值；‎‎2、‎A的‎与数‎?k‎??‎1?‎k‎?5‎01‎??‎1‎40‎最‎接近‎的特‎征值‎?i‎k（‎k=‎1，‎2，‎..‎.，‎39‎）；‎‎3、‎A的‎（谱‎范数‎）条‎件数‎c n‎d(‎A)‎2和‎行列‎式d‎e t‎A。

‎‎二、‎算法‎设计‎方案‎题‎目所‎给的‎矩阵‎阶数‎过大‎，必‎须经‎过去‎零压‎缩后‎进行‎存储‎和运‎算，‎本算‎法中‎压缩‎后的‎矩阵‎A1‎如下‎所示‎。

‎?0‎?0‎?A‎1?‎?a‎1‎??‎b?‎?c‎0‎b a‎2b‎c‎c b‎b c‎.‎..‎..‎..‎..‎..‎.‎c b‎b c‎c‎b a‎50‎0b‎0‎a ‎3.‎..‎a4‎99‎c‎?‎b?‎?a‎50‎1?‎?‎0?‎0?‎?‎由矩‎阵A‎的特‎征值‎满足‎的条‎件可‎知‎?1‎与?‎50‎1之‎间必‎有一‎个最‎大，‎则采‎用幂‎法求‎出的‎一‎个特‎征值‎必为‎其中‎的一‎个：‎当‎所求‎得的‎特征‎值为‎正数‎，则‎为?‎50‎1；‎否则‎为?‎1。

‎在求‎得?‎1与‎?‎50‎1其‎中的‎一个‎后，‎采用‎带位‎移的‎幂法‎则可‎求出‎它们‎中的‎另一‎个，‎且位‎移量‎即为‎先求‎出的‎特‎征值‎的值‎。

用‎反幂‎法求‎得的‎特征‎值必‎为?‎s。

‎由条‎件数‎的性‎质可‎得，‎c n‎d(‎A)‎2为‎模最‎大的‎特征‎值与‎模最‎小的‎特征‎值之‎比的‎模，‎因此‎，求‎出?‎1，‎?5‎01‎和?‎s的‎值后‎，则‎可以‎求得‎c n‎d(‎A)‎2。

北航计软实验报告实验三

北航计软实验报告实验三
3）针对所有的元素重复以上的步骤,除了最后一个。
4）持续每次对越来越少的元素重复上面的步骤,直到没有任何一对数字需要比较。
2.快速排序
1）设置两个变量i、j,排序开始的时候: i=0, j=N-1；
2）以第一个数组元素作为关键数据,赋值给key,即key=A[0]；
3）从j开始向前搜索,即由后开始向前搜索(j--),找到第一个小于key的值A[j],将A[j]和A[i]互换；
实验报告
实验名称冒泡排序和快速排序
班级
学号
姓名
成绩
实验概述:
【实验目的及要求】
通过编程程序达到熟悉并掌握教材中所介绍的几种排序方法。
【实验原理】
1.冒泡排序
1）比较相邻的元素。如果第一个比第二个大,就交换他们两个。
2）对每一对相邻元素作同样的工作,从开始第一对到结尾的最后一对。在这一点,最后的元素应该会是最大的数。
b[i]=num;
}
printf("随机数产生完成\n");
for (i=0;i<N;i++)
printf("%d ",a[i]);
printf("\n\n");
system("pause");
printf("\n\n");
com=0;
mov=0;
for (i=0;i<N;i++)
{
for (j=0;j<N-1;j++)
{
l++;
com++;
}
b[r]=b[l];
mov++;

数值分析上机实验报告

数值分析上机实验报告

数值分析上机实验报告摘要：本报告是对数值分析课程上机实验的总结和分析，涵盖了多种算法和数据处理方法，通过对实验结果的分析，探究了数值计算的一般过程和计算的稳定性。

1. 引言数值计算是数学的一个重要分支，广泛应用于物理、金融、工程等领域。

本次实验是对数值分析课程知识的实际应用，通过上机实现算法，探究数值计算的可靠性和误差分析。

2. 实验方法本次实验中，我们实现了多种算法，包括：（1）牛顿迭代法求方程的根；（2）高斯消元法求线性方程组的解；（3）最小二乘法拟合数据点；（4）拉格朗日插值法估计函数值；（5）梯形公式和辛普森公式求积分近似值。

对于每个算法，我们都进行了多组数值和不同参数的实验，并记录了相关数据和误差。

在实验过程中，我们着重考虑了算法的可靠性和计算的稳定性。

3. 实验结果与分析在实验中，我们得到了大量的实验数据和误差分析，通过对数据的展示和分析，我们得到了以下结论：（1）牛顿迭代法求解非线性方程的根能够对算法的初始值和迭代次数进行适当的调整，从而达到更高的稳定性和可靠性。

（2）高斯消元法求解线性方程组的解需要注意到矩阵的奇异性和精度的影响，从而保证计算的准确性。

（3）最小二乘法拟合数据点需要考虑到拟合的函数形式和数据的误差范围，采取适当的数据预处理和拟合函数的选择能够提高计算的准确性。

（4）拉格朗日插值法估计函数值需要考虑到插值点的选择和插值函数的阶数，防止出现龙格现象和插值误差过大的情况。

（5）梯形公式和辛普森公式求积分近似值需要考虑到采样密度和拟合函数的选择，从而保证计算的稳定性和收敛速度。

4. 结论通过本次实验的分析和总结，我们得到了深入的认识和理解数值计算的一般过程和算法的稳定性和可靠性，对于以后的数值计算应用也提供了一定的指导和参考。

北航数值分析报告大作业第三题(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 SY1415215DIMENSION X(11),Y(21),T(6),U(6),Z(6,6),UX(11,21),TY(11,21),FXY(11,21),C(6,6) DIMENSION X1(8),Y1(5),FXY1(8,5),PXY1(8,5),UX1(8,5),TY1(8,5)REAL(8) X,Y,T,U,Z,FXY,UX,TY,C,E,X1,Y1,FXY1,PXY1,UX1,TY1OPEN (1,FILE='第三题计算结果.TXT')DO I=1,11X(I)=0.08*(I-1)ENDDODO I=1,21Y(I)=0.5+0.05*(I-1)ENDDO!*****求解非线性方程组，得到z=f（t，u）的函数*******DO I=1,11DO J=1,21CALL DISNEWTON_NONLINEAR(X(I),Y(J),UX(I,J),TY(I,J)) ENDDOENDDO!*************分片二次插值得到z=f（x，y）***********DO I=1,11DO J=1,21CALL INTERPOLATION(UX(I,J),TY(I,J),FXY(I,J))ENDDOENDDOWRITE (1,"('数表(x,y,f(x,y)):')")WRITE (1,"(3X,'X',7X,'Y',10X,'F(X,Y)')")DO I=1,11DO J=1,21WRITE(1,'(1X,F5.2,2X,F5.3,2X,E20.13)') X(I),Y(J),FXY(I,J) ENDDOWRITE (1,"('')")ENDDO!***********最小二乘拟合得到P（x，y）**************N=11M=21WRITE (1,'(" ","K和σ分别为：")')DO K=1,20CALL LSFITTING(X,Y,FXY,C,N,M,K,K,E)WRITE (1,'(I3,2X,E20.13)') K-1,EIF(E<10E-7) EXITENDDOWRITE(1,'(" ")')WRITE(1,'("系数矩阵Crs（按行）为：")')DO I=1,KDO J=1,KWRITE (1,'(E20.13,2X,\)') C(I,J)ENDDOWRITE (1,"('')")WRITE (*,"('')")ENDDODO I=1,8X1(I)=0.1*IENDDODO J=1,5Y1(J)=0.5+0.2*JENDDODO I=1,8DO J=1,5CALL DISNEWTON_NONLINEAR(X1(I),Y1(J),UX1(I,J),TY1(I,J))ENDDOENDDODO I=1,8DO J=1,5CALL INTERPOLATION(UX1(I,J),TY1(I,J),FXY1(I,J))ENDDOENDDOPXY1=0DO I=1,8DO J=1,5DO II=1,KDO JJ=1,KPXY1(I,J)=PXY1(I,J)+C(II,JJ)*(X1(I)**(II-1))*(Y1(J)**(JJ-1)) ENDDOENDDOENDDOENDDOWRITE(1,'(" ")')WRITE(1,'("数表（x,y,f(x,y),p(x,y)):")')WRITE(1,"(2X,'X',6X,'Y',12X,'F(X,Y)',14X,'P(X,Y)')")DO I=1,8DO J=1,5WRITE(1,'(F5.3,2X,F5.3,2X,E20.13,2X,E20.13)') X1(I),Y1(J),FXY1(I,J),PXY1(I,J) ENDDOWRITE (1,"('')")ENDDOCLOSE (1)END!***********用离散牛顿法求解非线性方程组****************SUBROUTINE DISNEWTON_NONLINEAR(X1,Y1,U,T)PARAMETER (N=4)REAL EPS !EPS为迭代精度，M为最大迭代次数DIMENSION X(N),H(N),Y(N),JA(N,N),E(N),XK(N)REAL(8) JA,X,H,Y,E,XK,U,T,V,W,X1,Y1,E1,E2F1(T,U,V,W)=0.5*COS(T)+U+V+W-X1-2.67F2(T,U,V,W)=T+0.5*SIN(U)+V+W-Y1-1.07F3(T,U,V,W)=0.5*T+U+COS(V)+W-X1-3.74F4(T,U,V,W)=T+0.5*U+V+SIN(W)-Y1-0.79EPS=10E-12M=100X=1.0DO K=1,MH=1!计算Y=F（x）Y(1)=F1(X(1),X(2),X(3),X(4))Y(2)=F2(X(1),X(2),X(3),X(4))Y(3)=F3(X(1),X(2),X(3),X(4))Y(4)=F4(X(1),X(2),X(3),X(4))!计算JA（N,N）E=0.0DO I=1,NDO J=1,NDO JJ=1,NIF(JJ==J) THENE(JJ)=X(JJ)+H(JJ)ELSEE(JJ)=X(JJ)ENDIFENDDOIF(I==1) THENJA(I,J)=(F1(E(1),E(2),E(3),E(4))-F1(X(1),X(2),X(3),X(4)))/H(J) ELSEIF(I==2) THENJA(I,J)=(F2(E(1),E(2),E(3),E(4))-F2(X(1),X(2),X(3),X(4)))/H(J) ELSEIF(I==3) THENJA(I,J)=(F3(E(1),E(2),E(3),E(4))-F3(X(1),X(2),X(3),X(4)))/H(J) ELSEIF(I==4) THENJA(I,J)=(F4(E(1),E(2),E(3),E(4))-F4(X(1),X(2),X(3),X(4)))/H(J) ENDIFENDDOENDDO!求解线性方程组CALL GAUSS(JA,XK,-Y,N)!判断精度CALL NORM(XK,N,E1)CALL NORM(X,N,E2)IF(E1/E2<=EPS) THENT=X(1)U=X(2)EXITELSEDO I=1,NX(I)=X(I)+XK(I)ENDDOENDIFENDDORETURNEND!**********列主元高斯消去法求解线性方程组********* SUBROUTINE GAUSS(A,X,B,N)DIMENSION A(N,N),B(N),X(N),T(N,N),TB(N)REAL M(N,N)REAL(8) A,B,X,T!消元过程DO K=1,N-1TA=A(K,K)TL=KDO L=K+1,NIF ((A(L,K)>TA).OR.(A(L,K)==TA)) THENTA=A(L,K)TL=LDO J=K,NT(K,J)=A(K,J)A(K,J)=A(TL,J)A(TL,J)=T(K,J)ENDDOTB(K)=B(K)B(K)=B(TL)B(TL)=TB(K)ENDIFENDDODO I=K+1,NM(I,K)=A(I,K)/A(K,K)A(I,K)=0DO J=K+1,NA(I,J)=A(I,J)-M(I,K)*A(K,J)ENDDOB(I)=B(I)-M(I,K)*B(K)ENDDOENDDO!回代过程X(N)=B(N)/A(N,N)DO K=N-1,1,-1S=0.0DO J=K+1,NS=S+A(K,J)*X(J)ENDDOX(K)=(B(K)-S)/A(K,K)ENDDORETURNEND!***********求向量的无穷数************ SUBROUTINE NORM(X,N,A)DIMENSION X(N)REAL(8) X,AA=ABS(X(1))DO I=2,NIF(ABS(X(I))>ABS(X(I-1))) THENA=ABS(X(I))ENDIFENDDORETURNEND!**************分片二次代数插值************** SUBROUTINE INTERPOLATION(U,V,W) PARAMETER (N=6,M=6)DIMENSION X(N),Y(M),Z(M,N),LK(3),LR(3)REAL(8) X,Y,Z,H,TREAL(8) U,V,W,LK,LR !U,V分别为插值点处的坐标，W为插值结果INTEGER R!**********************数据赋值********************** DATA Y/0.0,0.2,0.4,0.6,0.8,1.0/DATA X/0.0,0.4,0.8,1.2,1.6,2.0/DATA Z/-0.5,-0.42,-0.18,0.22,0.78,1.5,&&-0.34,-0.5,-0.5,-0.34,-0.02,0.46,&&0.14,-0.26,-0.5,-0.58,-0.5,-0.26,&&0.94,0.3,-0.18,-0.5,-0.66,-0.66,&&2.06,1.18,0.46,-0.1,-0.5,-0.74,&&3.5,2.38,1.42,0.62,-0.02,-0.5/H=0.4T=0.2!******************计算K,R*************************IF(U<=X(2)+H/2) THENK=2ELSEIF(U>X(N-1)-H/2) THENK=N-1ELSEDO I=3,N-2IF((U>X(I)-H/2).AND.(U<=X(I)+H/2)) THENK=IENDIFENDDOENDIFIF(V<=Y(2)+T/2) THENR=2ELSEIF(V>Y(M-1)-T/2) THENR=M-1ELSEDO J=3,M-2IF((V>Y(J)-T/2).AND.(V<=Y(J)+T/2)) THENR=JENDIFENDDOENDIFI=KJ=RLK(1)=(U-X(I))*(U-X(I+1))/(X(I-1)-X(I))/(X(I-1)-X(I+1))LK(2)=(U-X(I-1))*(U-X(I+1))/(X(I)-X(I-1))/(X(I)-X(I+1)) LK(3)=(U-X(I))*(U-X(I-1))/(X(I+1)-X(I))/(X(I+1)-X(I-1)) LR(1)=(V-Y(J))*(V-Y(J+1))/(Y(J-1)-Y(J))/(Y(J-1)-Y(J+1)) LR(2)=(V-Y(J-1))*(V-Y(J+1))/(Y(J)-Y(J-1))/(Y(J)-Y(J+1)) LR(3)=(V-Y(J))*(V-Y(J-1))/(Y(J+1)-Y(J))/(Y(J+1)-Y(J-1))W=0DO K=1,3DO R=1,3W=W+LK(K)*LR(R)*Z(J+R-2,I+K-2)ENDDOENDDORETURNEND!*******************最小二乘拟合子函数************** SUBROUTINE LSFITTING(X,Y,Z,A,N,M,P,Q,DT1)INTEGER P,QDIMENSION X(N),Y(M),Z(N,M),A(P,Q)DIMENSION APX(20),APY(20),BX(20),BY(20),U(20,20),V(20,M) DIMENSION T(20),T1(20),T2(20)REAL(8) X,Y,Z,A,DT1DO I=1,PDO J=1,QA(I,J)=0.0ENDDOENDDOIF(P>N) P=NIF(P>20) P=20IF(Q>M) Q=MIF(Q>20) Q=20XX=0YY=0D1=NAPX(1)=0.0DO I=1,NAPX(1)=APX(1)+X(I)ENDDOAPX(1)=APX(1)/D1DO J=1,MV(1,J)=0.0DO I=1,NV(1,J)=V(1,J)+Z(I,J)ENDDOV(1,J)=V(1,J)/D1ENDDOIF(P>1) THEND2=0.0APX(2)=0.0DO I=1,NG=X(I)-APX(1)D2=D2+G*GAPX(2)=APX(2)+(X(I)-XX)*G*G ENDDOAPX(2)=APX(2)/D2BX(2)=D2/D1DO J=1,MV(2,J)=0.0DO I=1,NG=X(I)-APX(1)V(2,J)=V(2,J)+Z(I,J)*GENDDOV(2,J)=V(2,J)/D2ENDDOD1=D2ENDIFDO K=3,PD2=0.0APX(K)=0.0DO J=1,MV(K,J)=0.0ENDDODO I=1,NG1=1.0G2=X(I)-APX(1)DO J=3,KG=(X(I)-APX(J-1))*G2-BX(J-1)*G1 G1=G2G2=GENDDOD2=D2+G*GAPX(K)=APX(K)+X(I)*G*GDO J=1,MV(K,J)=V(K,J)+Z(I,J)*GENDDOENDDODO J=1,MV(K,J)=V(K,J)/D2ENDDOAPX(K)=APX(K)/D2BX(K)=D2/D1D1=D2ENDDOD1=MAPY(1)=0.0DO I=1,MAPY(1)=APY(1)+Y(I)ENDDOAPY(1)=APY(1)/D1DO J=1,PU(J,1)=0.0DO I=1,MU(J,1)=U(J,1)+V(J,I) ENDDOU(J,1)=U(J,1)/D1ENDDOIF(Q>1)THEND2=0.0APY(2)=0.0DO I=1,MG=Y(I)-APY(1)D2=D2+G*GAPY(2)=APY(2)+(Y(I))*G*G ENDDOAPY(2)=APY(2)/D2BY(2)=D2/D1DO J=1,PU(J,2)=0.0DO I=1,MG=Y(I)-APY(1)U(J,2)=U(J,2)+V(J,I)*GENDDOU(J,2)=U(J,2)/D2ENDDOD1=D2ENDIFDO K=3,QD2=0.0APY(K)=0.0DO J=1,PU(J,K)=0.0ENDDODO I=1,MG1=1.0G2=Y(I)-APY(1)DO J=3,KG=(Y(I)-APY(J-1))*G2-BY(J-1)*G1 G1=G2G2=GENDDOD2=D2+G*GAPY(K)=APY(K)+Y(I)*G*GDO J=1,PU(J,K)=U(J,K)+V(J,I)*GENDDOENDDODO J=1,PU(J,K)=U(J,K)/D2ENDDOAPY(K)=APY(K)/D2BY(K)=D2/D1D1=D2ENDDOV(1,1)=1.0V(2,1)=-APY(1)V(2,2)=1.0DO I=1,PDO J=1,QA(I,J)=0.0ENDDOENDDODO I=3,QV(I,I)=V(I-1,I-1)V(I,I-1)=-APY(I-1)*V(I-1,I-1)+V(I-1,I-2)IF(I>=4) THENDO K=I-2,2,-1V(I,K)=-APY(I-1)*V(I-1,K)+V(I-1,K-1)-BY(I-1)*V(I-2,K) ENDDOENDIFV(I,1)=-APY(I-1)*V(I-1,1)-BY(I-1)*V(I-2,1)ENDDODO I=1,PIF(I==1) THENT(1)=1.0T1(1)=1.0ELSEIF(I==2) THENT(1)=-APX(1)T(2)=1.0T2(1)=T(1)T2(2)=T(2)ELSET(I)=T2(I-1)T(I-1)=-APX(I-1)*T2(I-1)+T2(I-2)IF(I>=4) THENDO K=I-2,2,-1T(K)=-APX(I-1)*T2(K)+T2(K-1)-BX(I-1)*T1(K)ENDDOENDIFT(1)=-APX(I-1)*T2(1)-BX(I-1)*T1(1)T2(I)=T(I)DO K=I-1,1,-1T1(K)=T2(K)T2(K)=T(K)ENDDOENDIFDO J=1,QDO K=I,1,-1DO L=J,1,-1A(K,L)=A(K,L)+U(I,J)*T(K)*V(J,L)ENDDOENDDOENDDOENDDODT1=0.0DO I=1,NX1=X(I)DO J=1,MY1=Y(J)X2=1.0DD=0.0DO K=1,PG=A(K,Q)DO KK=Q-1,1,-1G=G*Y1+A(K,KK)ENDDOG=G*X2DD=DD+GX2=X2*X1ENDDODT=DD-Z(I,J)DT1=DT1+DT*DTENDDOENDDORETURNEND三、计算结果数表(x,y,f(x,y)):X Y UX TY F(X,Y) 0.00 0.500 1.345 0.243 0.17E+000.00 0.550 1.322 0.269 0.66E+000.00 0.600 1.299 0.295 0.35E+000.00 0.650 1.277 0.322 0.94E+000.00 0.700 1.255 0.350 0.30E-020.00 0.750 1.235 0.377 -0.87E-010.00 0.800 1.215 0.406 -0.58E+000.00 0.850 1.196 0.434 -0.72E+000.00 0.900 1.177 0.463 -0.54E+000.00 0.950 1.159 0.492 -0.86E+000.00 1.050 1.125 0.550 -0.74E+00 0.00 1.100 1.109 0.580 -0.06E+00 0.00 1.150 1.093 0.609 -0.00E+00 0.00 1.200 1.079 0.639 -0.18E+00 0.00 1.250 1.064 0.669 -0.52E+00 0.00 1.300 1.050 0.699 -0.19E+00 0.00 1.350 1.037 0.729 -0.48E+00 0.00 1.400 1.024 0.759 -0.68E+00 0.00 1.450 1.011 0.790 -0.52E+00 0.00 1.500 1.000 0.820 -0.29E+000.08 0.500 1.415 0.228 0.67E+00 0.08 0.550 1.391 0.253 0.08E+00 0.08 0.600 1.368 0.279 0.02E+00 0.08 0.650 1.346 0.306 0.47E+00 0.08 0.700 1.325 0.333 0.57E+00 0.08 0.750 1.304 0.360 0.48E-01 0.08 0.800 1.284 0.388 -0.73E-01 0.08 0.850 1.265 0.416 -0.16E+00 0.08 0.900 1.246 0.444 -0.29E+00 0.08 0.950 1.229 0.473 -0.36E+00 0.08 1.000 1.211 0.502 -0.08E+00 0.08 1.050 1.194 0.531 -0.29E+00 0.08 1.100 1.178 0.560 -0.78E+00 0.08 1.150 1.163 0.589 -0.93E+00 0.08 1.200 1.148 0.619 -0.44E+00 0.08 1.250 1.133 0.649 -0.92E+00 0.08 1.300 1.119 0.679 -0.71E+000.08 1.400 1.093 0.739 -0.37E+00 0.08 1.450 1.080 0.769 -0.83E+00 0.08 1.500 1.068 0.799 -0.92E+000.16 0.500 1.483 0.214 0.31E+00 0.16 0.550 1.460 0.239 0.64E+00 0.16 0.600 1.437 0.264 0.91E+00 0.16 0.650 1.414 0.290 0.06E+00 0.16 0.700 1.393 0.316 0.70E+00 0.16 0.750 1.372 0.343 0.59E+00 0.16 0.800 1.352 0.370 0.12E+00 0.16 0.850 1.333 0.398 0.77E-02 0.16 0.900 1.315 0.426 -0.83E-01 0.16 0.950 1.297 0.454 -0.58E+00 0.16 1.000 1.279 0.483 -0.20E+00 0.16 1.050 1.262 0.512 -0.11E+00 0.16 1.100 1.246 0.541 -0.74E+00 0.16 1.150 1.231 0.570 -0.09E+00 0.16 1.200 1.216 0.600 -0.59E+00 0.16 1.250 1.201 0.629 -0.66E+00 0.16 1.300 1.187 0.659 -0.71E+00 0.16 1.350 1.174 0.689 -0.32E+00 0.16 1.400 1.161 0.718 -0.56E+00 0.16 1.450 1.148 0.748 -0.31E+00 0.16 1.500 1.136 0.778 -0.75E+000.24 0.500 1.551 0.201 0.66E+01 0.24 0.550 1.527 0.225 0.03E+000.24 0.650 1.482 0.275 0.64E+00 0.24 0.700 1.460 0.301 0.47E+00 0.24 0.750 1.439 0.327 0.34E+00 0.24 0.800 1.419 0.354 0.24E+00 0.24 0.850 1.400 0.381 0.69E+00 0.24 0.900 1.381 0.409 0.04E-01 0.24 0.950 1.363 0.437 -0.42E-01 0.24 1.000 1.346 0.465 -0.06E+00 0.24 1.050 1.329 0.494 -0.59E+00 0.24 1.100 1.313 0.523 -0.83E+00 0.24 1.150 1.297 0.552 -0.15E+00 0.24 1.200 1.282 0.581 -0.19E+00 0.24 1.250 1.267 0.610 -0.84E+00 0.24 1.300 1.253 0.640 -0.66E+00 0.24 1.350 1.240 0.669 -0.30E+00 0.24 1.400 1.227 0.699 -0.86E+00 0.24 1.450 1.214 0.729 -0.84E+00 0.24 1.500 1.202 0.759 -0.77E+000.32 0.500 1.617 0.188 0.28E+01 0.32 0.550 1.593 0.212 0.49E+01 0.32 0.600 1.570 0.236 0.68E+00 0.32 0.650 1.547 0.261 0.75E+00 0.32 0.700 1.526 0.286 0.60E+00 0.32 0.750 1.505 0.312 0.77E+00 0.32 0.800 1.485 0.339 0.05E+00 0.32 0.850 1.466 0.365 0.99E+00 0.32 0.900 1.447 0.393 0.27E+000.32 1.000 1.411 0.448 -0.01E-02 0.32 1.050 1.395 0.477 -0.41E-01 0.32 1.100 1.378 0.505 -0.18E+00 0.32 1.150 1.363 0.534 -0.25E+00 0.32 1.200 1.347 0.563 -0.29E+00 0.32 1.250 1.333 0.592 -0.90E+00 0.32 1.300 1.319 0.621 -0.00E+00 0.32 1.350 1.305 0.650 -0.40E+00 0.32 1.400 1.292 0.680 -0.54E+00 0.32 1.450 1.279 0.710 -0.79E+00 0.32 1.500 1.267 0.739 -0.91E+000.40 0.500 1.681 0.177 0.91E+01 0.40 0.550 1.658 0.199 0.00E+01 0.40 0.600 1.634 0.223 0.83E+01 0.40 0.650 1.612 0.247 0.02E+01 0.40 0.700 1.591 0.272 0.94E+00 0.40 0.750 1.570 0.298 0.49E+00 0.40 0.800 1.550 0.324 0.94E+00 0.40 0.850 1.530 0.350 0.40E+00 0.40 0.900 1.512 0.377 0.33E+00 0.40 0.950 1.493 0.405 0.99E+00 0.40 1.000 1.476 0.432 0.68E+00 0.40 1.050 1.459 0.460 0.08E-01 0.40 1.100 1.443 0.488 -0.84E-01 0.40 1.150 1.427 0.517 -0.98E+00 0.40 1.200 1.412 0.545 -0.27E+00 0.40 1.250 1.397 0.574 -0.06E+000.40 1.350 1.369 0.632 -0.66E+00 0.40 1.400 1.356 0.662 -0.37E+00 0.40 1.450 1.343 0.691 -0.43E+00 0.40 1.500 1.331 0.721 -0.12E+000.48 0.500 1.745 0.166 0.69E+01 0.48 0.550 1.721 0.188 0.02E+01 0.48 0.600 1.698 0.211 0.74E+01 0.48 0.650 1.676 0.235 0.40E+01 0.48 0.700 1.654 0.259 0.23E+01 0.48 0.750 1.634 0.284 0.56E+00 0.48 0.800 1.613 0.310 0.28E+00 0.48 0.850 1.594 0.336 0.49E+00 0.48 0.900 1.575 0.363 0.31E+00 0.48 0.950 1.557 0.390 0.66E+00 0.48 1.000 1.539 0.417 0.30E+00 0.48 1.050 1.522 0.444 0.34E+00 0.48 1.100 1.506 0.472 0.07E-01 0.48 1.150 1.490 0.500 -0.62E-01 0.48 1.200 1.475 0.529 -0.45E+00 0.48 1.250 1.460 0.557 -0.86E+00 0.48 1.300 1.446 0.586 -0.39E+00 0.48 1.350 1.432 0.615 -0.22E+00 0.48 1.400 1.419 0.644 -0.67E+00 0.48 1.450 1.406 0.674 -0.55E+00 0.48 1.500 1.394 0.703 -0.14E+000.56 0.500 1.808 0.156 0.48E+010.56 0.600 1.761 0.200 0.10E+01 0.56 0.650 1.739 0.223 0.68E+01 0.56 0.700 1.717 0.247 0.94E+01 0.56 0.750 1.696 0.272 0.33E+01 0.56 0.800 1.676 0.297 0.11E+00 0.56 0.850 1.657 0.323 0.63E+00 0.56 0.900 1.638 0.349 0.97E+00 0.56 0.950 1.620 0.375 0.52E+00 0.56 1.000 1.602 0.402 0.56E+00 0.56 1.050 1.585 0.429 0.47E+00 0.56 1.100 1.568 0.457 0.20E+00 0.56 1.150 1.552 0.485 0.13E+00 0.56 1.200 1.537 0.513 0.09E-01 0.56 1.250 1.522 0.541 -0.47E-01 0.56 1.300 1.508 0.570 -0.99E+00 0.56 1.350 1.494 0.599 -0.82E+00 0.56 1.400 1.481 0.627 -0.26E+00 0.56 1.450 1.468 0.657 -0.71E+00 0.56 1.500 1.455 0.686 -0.98E+000.64 0.500 1.870 0.147 0.74E+01 0.64 0.550 1.846 0.168 0.10E+01 0.64 0.600 1.823 0.190 0.54E+01 0.64 0.650 1.801 0.213 0.42E+01 0.64 0.700 1.779 0.236 0.56E+01 0.64 0.750 1.758 0.260 0.03E+01 0.64 0.800 1.738 0.285 0.42E+01 0.64 0.850 1.718 0.310 0.41E+010.64 0.950 1.681 0.362 0.36E+00 0.64 1.000 1.664 0.388 0.18E+00 0.64 1.050 1.646 0.415 0.28E+00 0.64 1.100 1.630 0.443 0.07E+00 0.64 1.150 1.614 0.470 0.66E+00 0.64 1.200 1.598 0.498 0.09E+00 0.64 1.250 1.584 0.526 0.50E-01 0.64 1.300 1.569 0.554 -0.88E-01 0.64 1.350 1.555 0.583 -0.76E+00 0.64 1.400 1.542 0.611 -0.66E+00 0.64 1.450 1.529 0.640 -0.33E+00 0.64 1.500 1.516 0.669 -0.56E+000.72 0.500 1.931 0.139 0.94E+01 0.72 0.550 1.907 0.159 0.84E+01 0.72 0.600 1.884 0.181 0.36E+01 0.72 0.650 1.862 0.203 0.40E+01 0.72 0.700 1.840 0.226 0.47E+01 0.72 0.750 1.819 0.249 0.56E+01 0.72 0.800 1.799 0.273 0.19E+01 0.72 0.850 1.779 0.298 0.37E+01 0.72 0.900 1.760 0.323 0.86E+01 0.72 0.950 1.742 0.349 0.76E+00 0.72 1.000 1.724 0.375 0.24E+00 0.72 1.050 1.707 0.402 0.55E+00 0.72 1.100 1.691 0.429 0.97E+00 0.72 1.150 1.675 0.456 0.27E+00 0.72 1.200 1.659 0.484 0.31E+000.72 1.250 1.644 0.511 0.51E+00 0.72 1.300 1.630 0.539 0.49E+00 0.72 1.350 1.616 0.568 0.72E-02 0.72 1.400 1.602 0.596 -0.69E-01 0.72 1.450 1.589 0.625 -0.67E+00 0.72 1.500 1.576 0.653 -0.20E+000.80 0.500 1.992 0.131 0.31E+01 0.80 0.550 1.968 0.151 0.44E+01 0.80 0.600 1.945 0.172 0.41E+01 0.80 0.650 1.922 0.193 0.45E+01 0.80 0.700 1.900 0.216 0.00E+01 0.80 0.750 1.879 0.239 0.10E+01 0.80 0.800 1.859 0.263 0.16E+01 0.80 0.850 1.840 0.287 0.52E+01 0.80 0.900 1.821 0.312 0.02E+01 0.80 0.950 1.802 0.337 0.38E+01 0.80 1.000 1.784 0.363 0.89E+01 0.80 1.050 1.767 0.389 0.28E+00 0.80 1.100 1.751 0.416 0.09E+00 0.80 1.150 1.734 0.443 0.23E+00 0.80 1.200 1.719 0.470 0.93E+00 0.80 1.250 1.704 0.498 0.15E+00 0.80 1.300 1.689 0.525 0.86E+00 0.80 1.350 1.675 0.553 0.64E+00 0.80 1.400 1.662 0.582 0.74E-01 0.80 1.450 1.649 0.610 -0.37E-01 0.80 1.500 1.636 0.638 -0.81E+00K和σ分别为：0 0.93E+031 0.61E+012 0.92E-023 0.53E-034 0.16E-055 0.77E-07系数矩阵Crs（按行）为：0.00E+01 -0.83E+01 0.56E+00 0.97E+00 -0.03E+00 0.70E-010.91E+01 -0.99E+00 -0.96E+01 0.17E+01 -0.66E+00 0.10E-010.77E+00 0.42E+01 -0.10E+00 -0.81E+00 0.81E+00 -0.62E-01-0.25E+00 -0.21E+00 0.97E+00 -0.18E+00 0.49E+00 -0.63E-010.34E+00 -0.56E+00 0.69E-01 0.51E+00 -0.77E-01 0.27E-01-0.94E-01 0.94E+00 -0.58E+00 0.69E-01 -0.50E-01 0.53E-02数表（x,y,f(x,y),p(x,y)):X Y F(X,Y) P(X,Y)0.100 0.700 0.58E+00 0.05E+000.100 1.100 -0.66E+00 -0.26E+00 0.100 1.300 -0.68E+00 -0.31E+00 0.100 1.500 -0.52E+00 -0.49E+000.200 0.700 0.54E+00 0.19E+00 0.200 0.900 -0.63E-01 -0.65E-01 0.200 1.100 -0.90E+00 -0.90E+00 0.200 1.300 -0.84E+00 -0.90E+00 0.200 1.500 -0.03E+00 -0.04E+000.300 0.700 0.82E+00 0.09E+00 0.300 0.900 0.48E+00 0.11E+00 0.300 1.100 -0.63E+00 -0.88E+00 0.300 1.300 -0.72E+00 -0.96E+00 0.300 1.500 -0.34E+00 -0.84E+000.400 0.700 0.79E+00 0.89E+00 0.400 0.900 0.56E+00 0.63E+00 0.400 1.100 -0.83E-01 -0.04E-01 0.400 1.300 -0.72E+00 -0.71E+00 0.400 1.500 -0.85E+00 -0.07E+000.500 0.700 0.56E+01 0.92E+01 0.500 0.900 0.51E+00 0.23E+00 0.500 1.100 0.59E+00 0.27E+00 0.500 1.300 -0.53E+00 -0.11E+00 0.500 1.500 -0.67E+00 -0.33E+000.600 0.900 0.14E+00 0.75E+00 0.600 1.100 0.19E+00 0.32E+00 0.600 1.300 -0.70E-01 -0.82E-01 0.600 1.500 -0.08E+00 -0.75E+000.700 0.700 0.89E+01 0.29E+01 0.700 0.900 0.91E+01 0.11E+01 0.700 1.100 0.60E+00 0.97E+00 0.700 1.300 0.22E-01 0.06E-01 0.700 1.500 -0.53E+00 -0.80E+000.800 0.700 0.09E+01 0.06E+01 0.800 0.900 0.32E+01 0.50E+01 0.800 1.100 0.03E+00 0.79E+00 0.800 1.300 0.25E+00 0.50E+00 0.800 1.500 -0.14E+00 -0.28E+00。

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

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

北航数值分析全部三次大作业第一次大作业是关于解线性方程组的数值方法。

我们被要求实现各种常用的线性方程组求解算法，例如高斯消元法、LU分解法和迭代法等。

我首先学习了这些算法的原理和实现方法，并借助Python编程语言编写了这些算法的代码。

在实验中，我们使用了不同规模和条件的线性方程组进行测试，并比较了不同算法的性能和精度。

通过这个作业，我深入了解了线性方程组求解的原理和方法，提高了我的编程和数值计算能力。

第二次大作业是关于数值积分的方法。

数值积分是数值分析中的重要内容，它可以用于计算曲线的长度、函数的面积以及求解微分方程等问题。

在这个作业中，我们需要实现不同的数值积分算法，例如矩形法、梯形法和辛普森法等。

我学习了这些算法的原理和实现方法，并使用Python编写了它们的代码。

在实验中，我们计算了不同函数的积分值，并对比了不同算法的精度和效率。

通过这个作业，我深入了解了数值积分的原理和方法，提高了我的编程和数学建模能力。

第三次大作业是关于常微分方程的数值解法。

常微分方程是数值分析中的核心内容之一，它可以用于描述众多物理、化学和生物现象。

在这个作业中，我们需要实现不同的常微分方程求解算法，例如欧拉法、龙格-库塔法和Adams法等。

我学习了这些算法的原理和实现方法，并使用Python编写了它们的代码。

在实验中，我们解决了一些具体的常微分方程问题，并比较了不同算法的精度和效率。

通过这个作业，我深入了解了常微分方程的原理和方法，提高了我的编程和问题求解能力。

总的来说，北航数值分析课程的三次大作业非常有挑战性，但也非常有意义。

通过这些作业，我在数值计算和编程方面得到了很大的提升，也更加深入地了解了数值分析的理论和方法。

虽然这些作业需要大量的时间和精力，但我相信这些努力将会对我未来的学习和工作产生积极的影响。

BUAA数值分析大作业三

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(学硕)

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

数值分析上机实验报告

数值分析上机实验报告

数值分析上机实验报告导言：本次上机实验主要是针对数值分析课程中的一些基本算法进行实验验证。

实验内容包括迭代法、插值法、数值积分和常微分方程的数值解等。

在实验过程中，我们将会使用MATLAB进行算法的实现，并对结果进行分析。

一、迭代法迭代法是解决函数零点、方程解等问题的常用方法。

我们将选择几个常见的函数进行迭代求根的实验。

（1）二分法二分法是一种简单而有效的迭代求根法。

通过函数在区间两个端点处的函数值异号来确定函数在区间内存在零点，并通过不断缩小区间来逼近零点。

（2）牛顿法牛顿法利用函数的一阶导数和二阶导数的信息来逼近零点。

通过不断迭代更新逼近值，可以较快地求得零点。

实验结果表明，对于简单的函数，这两种迭代法都具有很好的收敛性和稳定性。

但对于一些复杂的函数，可能会出现迭代失效或者收敛速度很慢的情况。

二、插值法插值法是在给定一些离散数据点的情况下，通过构造一个插值函数来逼近未知函数的值。

本实验我们将使用拉格朗日插值和牛顿插值两种方法进行实验。

（1）拉格朗日插值拉格朗日插值通过构造一个多项式函数来逼近未知函数的值。

该多项式经过离散数据点，并且是唯一的。

该方法简单易懂，但插值点越多，多项式次数越高，插值函数的精度也就越高。

（2）牛顿插值牛顿插值利用差商的概念，通过构造一个插值多项式来逼近未知函数的值。

与拉格朗日插值相比，牛顿插值的计算过程更加高效。

但同样要求插值点的选择要合理，否则可能出现插值函数不收敛的情况。

实验结果表明，这两种插值方法都能够很好地逼近未知函数的值。

插值点的选择对插值结果有很大的影响，过多或者过少的插值点都可能导致插值结果偏离真实函数的值。

三、数值积分数值积分是一种将定积分问题转化为数值求和的方法。

本实验我们将使用复合梯形求积法和复合辛普森求积法进行实验。

（1）复合梯形求积法复合梯形求积法将定积分区间等分为若干小区间，然后使用梯形公式对每个小区间进行近似求积，最后将结果相加得到整个定积分的近似值。

北航数分大作业三

北航数分大作业三

一、算法的设计方案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 。

北航数值分析大作业3

北航数值分析大作业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 integral implicit none contains!//////////复化梯形 subroutine trapezoid(m) implicit none integer :: i,j,k,mreal*8 :: x(m+1),u(m+1) real*8 :: sum,sum1,g,r,h real*8 :: e=1.0e-10h=2./m do i=1,m+1x(i)=-1.+(i-1)*h end dou=0.02 do k=1,50 do i=1,m+1 sum1=0.g=dexp(x(i)*4.)+(dexp(x(i)+4.)-dexp(-4.-x(i)))/(x(i)+4.) do j=2,m sum1=sum1+dexp(x(i)*x(j))*u(j) end do sum=h/2.*(dexp(x(i)*-1.)*u(1)+dexp(x(i)*1.)*u(m+1)+2*sum1) u(i)=g-sum end 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.各方法所得曲线（由于所取节点太多，且精度高，所以图中很难看出各曲线的区别。

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

北航数值分析大作业题目三
《数值分析》第三次大作业
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},

数值分析上机报告资料

数值分析上机报告资料

数值分析上机报告姓名：学号：专业：联系电话：本次数值分析上机实习采用Matlab数学软件。

Matlab是一种用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境。

在数值分析应用中可以直接调用Matlab软件中已有的函数，同时用户也可以将自己编写的实用导入到Matlab函数库中方便自己调用。

基于Matlab数学软件的各种实用性功能与优点，本次数值分析实习决定采用其作为分析计算工具。

1．语言简洁，编程效率高因为MATLAB定义了专门用于矩阵运算的运算符，使得矩阵运算就像列出算式执行标量运算一样简单，而且这些运算符本身就能执行向量和标量的多种运算。

利用这些运算符可使一般高级语言中的循环结构变成一个简单的MATLAB语句，再结合MATLAB丰富的库函数可使变得相当简短，几条语句即可代替数十行C语言或Fortran语言语句的功能。

2. 交互性好，使用方便在MATLAB的命令窗口中，输入一条命令，立即就能看到该命令的执行结果，体现了良好的交互性。

交互方式减少了编程和调试的工作量，给使用者带来了极大的方便。

因为不用像使用C语言和Fortran语言那样，首先编写源，然后对其进行编译、连接，待形成可执行文件后，方可运行得出结果。

3. 强大的绘图能力，便于数据可视化MATLAB不仅能绘制多种不同坐标系中的二维曲线，还能绘制三维曲面，体现了强大的绘图能力。

正是这种能力为数据的图形化表示(即数据可视化)提供了有力工具，使数据的展示更加形象生动，有利于揭示数据间的内在关系在新版本中也加入了对C、FORTRAN、c++、JA V A的支持，使用时可以直接调用，也可将编写的实用程序导入到matlab函数库中方便以后使用时调用。

本次编程所用的软件为MATLAB，通过这次作业，对它有了初步的认识，以及对数值分析的体会更为深刻，希望为以后的学习和工作奠定一定的基。

目录1 必做题一插值法 (4)1.1题目 (4)1.2 分析过程 (4)1.3 计算结果 (5)1.4 结果分析 (6)2 必做题二雅格比法迭代与高斯－赛德尔迭代 (6)2.1题目 (6)2.2分析过程 (6)2.3计算结果 (7)2.4 结果分析 (8)3 选做题一 (8)3.1题目三次样条插值 (8)3.2分析过程 (8)3.3计算结果 (9)3.4 结果分析 (9)附录 (10)附录一：必做题一插值法代码 (11)附录二：必做题二雅格比法迭代与高斯－赛德尔迭代代码 (12)附录三：选做题一三次样条插值代码 (14)1 必做题一 插值法1.1题目某过程涉及两变量x 和y, 拟分别用插值多项式和多项式拟合给出其对应规律的近似多项式，已知xi 与yi 之间的对应数据如下，xi=1,2,…,10yi = 34.6588 40.3719 14.6448 -14.2721 -13.3570 24.8234 75.2795 103.5743 97.4847 78.2392（1）请用次数分别为3，4，5，6的多项式拟合并给出最好近似结果f(x)。

北航数值分析A大作业3

北航数值分析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. 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

数值分析上机作业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 。

在拟合,0(,)kr s rsr s p x y cx y ==∑时，需要计算11()()T T T C B B B UG G G --=，令,T T W B B M G G ==，计算,W M 时，根据对称性只需要计算对角线元素和对角线以上元素即可，节省运算时间。

于是TWCM B UG E ==，用选主元的LU 分解法求解WF E =，再计算TTMC F =，这里,C F 只需要按行取元素进行运算即可，故不需要进行转置运算。

k 从1到9依次增加，计算σ的值，当710σ-≤时，得到达到精度的最小的k 。

2.打印输出结果f(x,y):4.46504018480e-001 3.24683262927e-001 2.10159686683e-001 1.03043608316e-0013.40189556266e-003 -8.87358136384e-0026.38015226510e-001 5.0661*******e-001 3.82176369277e-001 2.64863491154e-0011.54780200285e-001 5.199********e-0028.40081395765e-001 6.99764165673e-001 5.66061442351e-001 4.39171608118e-0013.19242138041e-001 2.06376192387e-0011.0515*******e+000 9.029********e-001 7.60580266859e-001 6.24715198145e-0014.95519756001e-001 3.73134042774e-0011.27124675148e+000 1.11500201815e+000 9.64607727215e-001 8.20347369475e-0016.82447678179e-001 5.51085208597e-0011.49832105248e+000 1.33499863207e+000 1.177********e+000 1.025********e+0008.78960023174e-001 7.39145108703e-0011.73189274038e+000 1.56203457721e+000 1.39721691821e+000 1.23780100674e+0001.0840*******e+000 9.36322772315e-0011.97122178640e+000 1.79532959950e+000 1.62406711323e+000 1.45783058271e+0001.29695464975e+000 1.14171810545e+000k=1, sigma=3.22090897363e+000k=2, sigma=4.65996003320e-003k=3, sigma=1.72117537926e-004k=4, sigma=3.30953430190e-006k=5, sigma=2.54137773513e-008C_rs:2.021********e+000 -3.66842591519e+000 7.09246688452e-001 8.48607444215e-001-4.158********e-001 6.74322022261e-0023.19192646165e+000 -7.41209812962e-001 -2.69690653026e+000 1.63095275395e+000-4.84601977266e-001 6.05908514261e-0022.56706343343e-001 1.58096413206e+000 -4.65701259544e-001 -7.89186706497e-0021.00853116187e-001 -2.0768*******e-002-2.68608304872e-001 -7.33963449843e-001 1.08429601112e+000 -8.156********e-0013.0728*******e-001 -4.68489486618e-0022.16521800059e-001 -1.73026852965e-001 -8.41324310602e-002 2.55736987891e-001-1.47683427939e-001 2.77711894906e-002-5.54328606191e-002 1.40518220408e-001 -1.30388672239e-001 3.44966421960e-0026.95959872154e-003 -3.30022941240e-003f(x*,y*):0.194720 -0.183037 -0.445498 -0.597567 -0.646460 0.405979 -0.022516 -0.338221 -0.544438 -0.647361 0.634777 0.158801 -0.207366 -0.465358 -0.6202710.878960 0.358651 -0.055253 -0.362680 -0.5675651.136611 0.574980 0.115992 -0.238568 -0.491434 1.406042 0.805941 0.304429 -0.095016 -0.393902 1.685784 1.049881 0.508294 0.066149 -0.276834 1.974575 1.305335 0.725992 0.243254 -0.141950p(x*,y*):0.194730 -0.183042 -0.445500 -0.597559 -0.646446 0.405990 -0.022521 -0.338224 -0.544430 -0.647348 0.634787 0.158796 -0.207369 -0.465350 -0.6202570.878970 0.358646 -0.055255 -0.362671 -0.5675511.136620 0.574976 0.115989 -0.238560 -0.491421 1.406051 0.805937 0.304426 -0.095009 -0.393890 1.685791 1.049878 0.508291 0.066156 -0.276822 1.974581 1.305332 0.725989 0.243261 -0.141939f x y插值的数值情况：对(,)(,)(,)i i i i f x y p x y 的数值情况（z 坐标为10-7）：3. 源程序 1) 主函数#include<stdio.h> #include<math.h>#define length_x 11//x 向量长度 #define length_y 21//y 向量长度 #define x_length 8//x*向量长度 #define y_length 5//y*向量长度extern double interpolation(double tt,double uu,double **table); //输入一点，输出一点插值extern double row_multi_sum(double **B,int i,int j,int length); //矩阵列相乘int main() {double tt=0; double uu=0;double table[6][6]={0}; double ff=0;double xx[length_x]={0};double yy[length_y]={0};double sum=0.0;double U[length_x][length_y]={0};double B[length_x][10]={0};double G[length_y][10]={0};double W[10][10]={0};double T[10][10]={0};double M[10][10]={0};double C[10][10]={0};double F[10][length_y]={0};double E[10][10]={0};double XX[10][10]={0};//double ff[length_x][length_y]={0}; //记录f阵double pp[length_x][length_y]={0}; //记录p阵int r=0;int s=0;int ii,jj,kk;int k_max=0;FILE *fd;fd = fopen("tmptmp.txt", "w");ii=0;jj=0;kk=0;table[0][0]=-0.5;table[0][1]=-0.34;table[0][2]=0.14;table[0][3]=0.94;table[0][4]=2.06;table[0][5]=3.5;table[1][0]=-0.42;table[1][1]=-0.5;table[1][2]=-0.26;table[1][3]=0.3;table[1][4]=1.18;table[1][5]=2.38;table[2][0]=-0.18;table[2][1]=-0.5;table[2][2]=-0.5;table[2][3]=-0.18;table[2][4]=0.46;table[2][5]=1.42;table[3][0]=0.22;table[3][1]=-0.34;table[3][2]=-0.58;table[3][3]=-0.5;table[3][4]=-0.1;table[3][5]=0.62;table[4][0]=0.78;table[4][1]=-0.02;table[4][2]=-0.5;table[4][3]=-0.66;table[4][4]=-0.5;table[4][5]=-0.02;table[5][0]=1.5;table[5][1]=0.46;table[5][2]=-0.26;table[5][3]=-0.66;table[5][4]=-0.74;table[5][5]=-0.5;for(ii=0;ii<length_x;ii++)xx[ii]=0.08*ii;for(ii=0;ii<length_y;ii++)yy[ii]=0.5+0.05*ii;for(ii=0;ii<length_x;ii++) //生成U阵{for(jj=0;jj<length_y;jj++){ite_equations(&tt,&uu,xx[ii],yy[jj]); //牛顿迭代解方程，输入一点（x,y）输出一点(t,u)U[ii][jj]=(double)interpolation(tt,uu,(double **)table); //输入一点，输出一点插值}}fprintf(fd,"\nf(x,y):\n");for(ii=0;ii<x_length;ii++){for(jj=0;jj<=y_length;jj++){fprintf(fd,"%.11e ",U[ii][jj]);}fprintf(fd,"\n");}for(k_max=1;k_max<=6;k_max++){for(kk=0;kk<(k_max+1);kk++) //生成B阵{for(ii=0;ii<length_x;ii++){B[ii][kk]=pow(xx[ii],kk);}}for(kk=0;kk<(k_max+1);kk++){for(ii=0;ii<length_y;ii++) //生成G阵{G[ii][kk]=pow(yy[ii],kk);}}for(ii=0;ii<=k_max;ii++) //BT * B{for(jj=0;jj<=k_max;jj++){W[ii][jj]=row_multi_sum((double **)B,ii,jj,length_x);}}for(ii=0;ii<=k_max;ii++) //GT * G{for(jj=0;jj<=k_max;jj++){M[ii][jj]=row_multi_sum((double **)G,ii,jj,length_y);}}for(ii=0;ii<=k_max;ii++) //F=B'*Ufor(jj=0;jj<length_y;jj++){F[ii][jj]=0;for(kk=0;kk<length_x;kk++)F[ii][jj]=F[ii][jj]+B[kk][ii]*U[kk][jj];}for(ii=0;ii<=k_max;ii++) //E=F*Gfor(jj=0;jj<=k_max;jj++){E[ii][jj]=0;for(kk=0;kk<length_y;kk++)E[ii][jj]=E[ii][jj]+F[ii][kk]*G[kk][jj];}Matrix_LU((double **)W,(double **)E,(double **)M,(double**)C,k_max,0);//(double **)XX,for(ii=0;ii<length_x;ii++)for(jj=0;jj<length_y;jj++){sum=0;for(r=0;r<=k_max;r++)for(s=0;s<=k_max;s++)sum=sum+C[r][s]*pow(xx[ii],r)*pow(yy[jj],s);pp[ii][jj]=sum;}sum=0;for(ii=0;ii<length_x;ii++)for(jj=0;jj<length_y;jj++){sum=sum+(pp[ii][jj]-U[ii][jj])*(pp[ii][jj]-U[ii][jj]);}fprintf(fd,"\nk=%d, sigma=%.11e\n",k_max,sum);if(sum<=1e-7){fprintf(fd,"\nC_rs:\n");for(r=0;r<=k_max;r++){for(s=0;s<=k_max;s++){fprintf(fd,"%.11e ",C[r][s]);}fprintf(fd,"\n");}break;}}for(ii=0;ii<x_length;ii++)xx[ii]=0.1*(ii+1);for(ii=0;ii<y_length;ii++)yy[ii]=0.5+0.2*(ii+1);for(ii=0;ii<x_length;ii++) //生成f*阵{for(jj=0;jj<y_length;jj++){ite_equations(&tt,&uu,xx[ii],yy[jj]); //牛顿迭代解方程，输入一点（x,y）输出一点(t,u)U[ii][jj]=(double)interpolation(tt,uu,(double **)table); //输入一点，输出一点插值}}fprintf(fd,"\nf(x*,y*):\n");for(ii=0;ii<x_length;ii++){for(jj=0;jj<y_length;jj++){fprintf(fd,"%-10f ",U[ii][jj]);}fprintf(fd,"\n");}for(ii=0;ii<x_length;ii++) //生成p*阵{for(jj=0;jj<y_length;jj++){sum=0;for(r=0;r<=k_max;r++)for(s=0;s<=k_max;s++)sum=sum+C[r][s]*pow(xx[ii],r)*pow(yy[jj],s);pp[ii][jj]=sum;}}fprintf(fd,"\np(x*,y*):\n");for(ii=0;ii<x_length;ii++){for(jj=0;jj<y_length;jj++){fprintf(fd,"%-10f ",pp[ii][jj]);}fprintf(fd,"\n");}fclose(fd);}2)矩阵2列相乘#include<stdio.h>#include<math.h>#define size 10#define TTT(a,b) *((double *)TTT+(a)*(size)+b) //L U 的宏定义，以方便操作二维数组double row_multi_sum(double **B,int i,int j,int length) //矩阵列相乘{double sum=0;double **TTT;int kk=0;TTT=B;for(kk=0;kk<length;kk++){sum=sum+(TTT(kk,i))*(TTT(kk,j));}return sum;}3)9点分片双二次插值#include<stdio.h>#include<math.h>#define size 6#define table(a,b) *((double *)table+(a)*(size)+b) //L U 的宏定义，以方便操作二维数组double interpolation(double tt,double uu,double **table) //输入一点，输出一点插值{int ii=0;int jj=0;int section_t=0;int section_u=0;double answer=0;double temp=0;double temp_u[3]={0};double temp_t[3]={0};double my_t[6]={0};double my_u[6]={0};for(ii=1;ii<6;ii++){my_u[ii]=my_u[ii-1]+0.4;my_t[ii]=my_t[ii-1]+0.2;}//判断输入点在哪个区域//判断u用哪段插值if ((uu>=0)&&(uu<0.6))section_u=1;else if((uu>=0.6)&&(uu<1.0))section_u=2;else if((uu>=1.0)&&(uu<1.4))section_u=3;else if((uu>=1.4)&&(uu<=2))section_u=4;else;//判断t用哪段插值if ((tt>=0)&&(tt<0.3))section_t=1;else if((tt>=0.3)&&(tt<0.5))section_t=2;else if((tt>=0.5)&&(tt<0.7))section_t=3;else if((tt>=0.7)&&(tt<=1.0))section_t=4;else;temp_u[0]=(uu-my_u[section_u])*(uu-my_u[section_u+1])/((my_u[section _u-1]-my_u[section_u])*(my_u[section_u-1]-my_u[section_u+1]));temp_u[1]=(uu-my_u[section_u-1])*(uu-my_u[section_u+1])/((my_u[secti on_u]-my_u[section_u-1])*(my_u[section_u]-my_u[section_u+1]));temp_u[2]=(uu-my_u[section_u-1])*(uu-my_u[section_u])/((my_u[section _u+1]-my_u[section_u-1])*(my_u[section_u+1]-my_u[section_u]));temp_t[0]=(tt-my_t[section_t])*(tt-my_t[section_t+1])/((my_t[section _t-1]-my_t[section_t])*(my_t[section_t-1]-my_t[section_t+1]));temp_t[1]=(tt-my_t[section_t-1])*(tt-my_t[section_t+1])/((my_t[secti on_t]-my_t[section_t-1])*(my_t[section_t]-my_t[section_t+1]));temp_t[2]=(tt-my_t[section_t-1])*(tt-my_t[section_t])/((my_t[section _t+1]-my_t[section_t-1])*(my_t[section_t+1]-my_t[section_t]));answer=0;for(ii=0;ii<3;ii++){answer=answer+temp_t[ii]*(temp_u[0]*table((section_t-1+ii),section_u -1)+temp_u[1]*table((section_t-1+ii),section_u)+temp_u[2]*table((sectio n_t-1+ii),section_u+1));}return answer;}4)Newton迭代求解非线性方程组#include<stdio.h>#include<math.h>void ite_equations(double *tt,double *uu,double xx,double yy) //牛顿迭代解非线性方程组{double t=1;double u=1;double v=1;double w=1;double d[4]={0}; //增量double x=xx;double y=yy;double F_Jacob[4][4]={0};double F[4]={0};double L[9][9]={0};double U[9][9]={0};double temp1=0;double temp2=0;int ii=0;int jj=0;int kk=0;for(ii=0;ii<=1e2;ii++){F_Jacob[0][0]=-0.5*sin(t);F_Jacob[0][1]=1;F_Jacob[0][2]=1;F_Jacob[0][3]=1;F_Jacob[1][0]=1;F_Jacob[1][1]=0.5*cos(u);F_Jacob[1][2]=1;F_Jacob[1][3]=1;F_Jacob[2][0]=0.5;F_Jacob[2][1]=1;F_Jacob[2][2]=-sin(v);F_Jacob[2][3]=1;F_Jacob[3][0]=1;F_Jacob[3][1]=0.5;F_Jacob[3][2]=1;F_Jacob[3][3]=cos(w);F[0]=0.5*cos(t)+u+v+w-x-2.67;F[1]=t+0.5*sin(u)+v+w-y-1.07;F[2]=0.5*t+u+cos(v)+w-x-3.74;F[3]=t+0.5*u+v+sin(w)-y-0.79;LU((double **)L,(double **)U,(double **)F_Jacob,(double **)F);huidai((double **)L,(double **)U,(double **)F,(double **)d);t=t+d[0];u=u+d[1];v=v+d[2];w=w+d[3];temp1=sqrt(t*t+u*u+v*v+w*w);temp2=sqrt(d[0]*d[0]+d[1]*d[1]+d[2]*d[2]+d[3]*d[3]);if((temp2/temp1)<1e-12){*tt=t;*uu=u;break;}}5)解线性方程组用到LU分解#include<stdio.h>#include<math.h>#define size 4//#define halfBW 3#define L(a,b) *((double *)L+(a)*(size)+b) //L U 的宏定义，以方便操作二维数组#define U(a,b) *((double *)U+(a)*(size)+b)#define A(a,b) *((double *)A+(a)*(size)+b) //A 的宏定义，以方便操作二维数组//#define B(a,b) *((double *)B+(a)*(size)+b) //B 的宏定义，以方便操作二维数组//#define X(a,b) *((double *)X+(a)*(size)+b) //A 的宏定义，以方便操作二维数组#define b(ii) *((double *)b+ii)void LU(double **LL,double **UU,double **AA,double **bb) // //int a //需要选主元，{double **L;double **U;double **A;double A_temp[size][size]={0};double **b;double y[size]={0};double x[size]={0};int kk,ii,jj,pp,qq;double s,s_temp=0;double sum_answer=0;double temp=0;int s_count=0;int n=size;////////////////L=LL;U=UU;A=AA;b=bb;/////////////kk=0;ii=0;jj=0;pp=0;qq=0;//////////////////////////////////for(ii=0;ii<n;ii++) //对A_temp赋值for(jj=0;jj<n;jj++)A_temp[ii][jj]=A(ii,jj);for(ii=0;ii<n;ii++)L(ii,ii)=1.0;for (kk=0;kk<n;kk++) //总循环求出L U{ sum_answer=0;s_count=kk;for(pp=kk;pp<n;pp++)//doolitte 选主元{if (kk==0)sum_answer=0;elsefor(qq=0;qq<kk;qq++)sum_answer=sum_answer+L(pp,qq)*U(qq,pp);s=A_temp[pp][kk]-sum_answer;if (s>s_temp){s_temp=s;s_count=pp;}else if (s<(-1.0)*s_temp){s_temp=(-1.0)*s;s_count=pp;}else;}for(pp=0;pp<n;pp++) //换位A_temp{temp=A_temp[s_count][pp];A_temp[s_count][pp]=A_temp[kk][pp];A_temp[kk][pp]=temp;}for(pp=0;pp<kk;pp++) //换位L{temp=L(s_count,pp);L(s_count,pp)=L(kk,pp);L(kk,pp)=temp;}temp=b(s_count);b(s_count)=b(kk);b(kk)=temp;for (ii=kk;ii<n;ii++) //先求一行U{if(kk==0)U(kk,ii)=A_temp[kk][ii];else{sum_answer=0;for(qq=0;qq<kk;qq++)sum_answer=sum_answer+L(kk,qq)*U(qq,ii);U(kk,ii)=A_temp[kk][ii]-sum_answer;}}for(ii=kk+1;ii<n;ii++) //求出一列L{sum_answer=0;for(qq=0;qq<kk;qq++)sum_answer=sum_answer+L(ii,qq)*U(qq,kk);L(ii,kk)=(A_temp[ii][kk]-sum_answer)/U(kk,kk);}}}6)回代求解方程组#include<stdio.h>#include<math.h>#define size 4#define L(a,b) *((double *)L+(a)*(size)+b) //L U 的宏定义，以方便操作二维数组#define U(a,b) *((double *)U+(a)*(size)+b)#define x(a) *((double *)x+a) //L U 的宏定义，以方便操作二维数组#define y(a) *((double *)y+a)#define b(a) *((double *)b+a)//#define xx(a) xx[a]//A*uu=yy => L*U*uu=yy =>L*xx=yy && U*uu=xxvoid huidai(double **LL,double **UU,double **bb,double **xx) //回代法求xx, LUx=b{double **L;double **U;double **b;double **x;double y[size]={0};double sum_answer=0;int ii=0;int jj=0;L=LL;U=UU;b=bb;x=xx;//迭代Ly=by[0]=-b(0);for(ii=1;ii<size;ii++){sum_answer=0;for(jj=1;jj<ii;jj++){sum_answer=sum_answer+L(ii,jj)*y[jj];}y[ii]=-b(ii)-sum_answer;}//迭代Ux=yx(size-1)=y[size-1]/U((size-1),(size-1));for(ii=size-2;ii>=0;ii--){sum_answer=0;for(jj=ii+1;jj<size;jj++){sum_answer=sum_answer+U(ii,jj)*x(jj);}x(ii)=(y[ii]-sum_answer)/U(ii,ii);}}7)求解C=（B’*B）-1*(B’*U*G)*(G’*G)#include<stdio.h>#include<math.h>#define size 10//#define halfBW 3#define L(a,b) L[a][b]//*((double *)L+(a)*(size)+b) //L U 的宏定义，以方便操作二维数组#define U(a,b) U[a][b]//*((double *)U+(a)*(size)+b)#define A(a,b) *((double *)A+(a)*(size)+b) //A 的宏定义，以方便操作二维数组#define B(a,b) *((double *)B+(a)*(size)+b) //B 的宏定义，以方便操作二维数组#define M(a,b) *((double *)M+(a)*(size)+b) //B 的宏定义，以方便操作二维数组#define C(a,b) *((double *)C+(a)*(size)+b) //B 的宏定义，以方便操作二维数组//#define X(a,b) *((double *)X+(a)*(size)+b) //A 的宏定义，以方便操作二维数组#define X(a,b) X[a][b]//#define b(ii) *((double *)b+ii)void Matrix_LU(double **AA,double **BB,double **MM,double **CC,intk_max,int mode){double **A;double **B;double **M;double **C;double X[size][size]={0};int n=k_max+1;double A_temp[size][size]={0}; double L[size][size]={0};double U[size][size]={0};double Y[size][size]={0};double **b;int kk,ii,jj,pp,qq;double s,s_temp=0;double sum_answer=0;double temp=0;int s_count=0;//////////////////L=LL;//U=UU;A=AA;B=BB;M=MM;C=CC;//X=XX;/////////////kk=0;ii=0;jj=0;pp=0;qq=0;////////////////////////////////// for(ii=0;ii<n;ii++) //对A_temp赋值for(jj=0;jj<n;jj++)A_temp[ii][jj]=A(ii,jj);for(ii=0;ii<n;ii++)L(ii,ii)=1.0;for (kk=0;kk<n;kk++) //总循环求出L U{ sum_answer=0;s_count=kk;for(pp=kk;pp<n;pp++)//doolitte 选主元{if (kk==0)sum_answer=0;elsefor(qq=0;qq<kk;qq++)sum_answer=sum_answer+L(pp,qq)*U(qq,pp);s=A_temp[pp][kk]-sum_answer;if (s>s_temp){s_temp=s;s_count=pp;}else if (s<(-1.0)*s_temp){s_temp=(-1.0)*s;s_count=pp;}else;}for(pp=0;pp<n;pp++) //换位A_temp{temp=A_temp[s_count][pp];A_temp[s_count][pp]=A_temp[kk][pp];A_temp[kk][pp]=temp;}for(pp=0;pp<kk;pp++) //换位L{temp=L(s_count,pp);L(s_count,pp)=L(kk,pp);L(kk,pp)=temp;}for(pp=0;pp<n;pp++) //换位B /////////////{temp=B(s_count,pp);B(s_count,pp)=B(kk,pp);B(kk,pp)=temp;}for (ii=kk;ii<n;ii++) //先求一行U{if(kk==0)U(kk,ii)=A_temp[kk][ii];else{sum_answer=0;for(qq=0;qq<kk;qq++)sum_answer=sum_answer+L(kk,qq)*U(qq,ii);U(kk,ii)=A_temp[kk][ii]-sum_answer;}}for(ii=kk+1;ii<n;ii++) //求出一列L{sum_answer=0;for(qq=0;qq<kk;qq++)sum_answer=sum_answer+L(ii,qq)*U(qq,kk);L(ii,kk)=(A_temp[ii][kk]-sum_answer)/U(kk,kk);}}//迭代Ly=bfor(kk=0;kk<n;kk++){Y[0][kk]=B(0,kk);for(ii=1;ii<n;ii++){sum_answer=0;for(jj=0;jj<ii;jj++){sum_answer=sum_answer+L(ii,jj)*Y[jj][kk];}Y[ii][kk]=B(ii,kk)-sum_answer;}//迭代Ux=yX((n-1),kk)=Y[n-1][kk]/U((n-1),(n-1));for(ii=n-2;ii>=0;ii--){sum_answer=0;for(jj=ii+1;jj<n;jj++){sum_answer=sum_answer+U(ii,jj)*X(jj,kk);}X(ii,kk)=(Y[ii][kk]-sum_answer)/U(ii,ii);}}////////////////////////////////////////////////////////////// // 再操作，（G'G）C'=X'//////////////////////////////////for(ii=0;ii<n;ii++) //对A_temp赋值for(jj=0;jj<n;jj++)A_temp[ii][jj]=M(ii,jj);for(ii=0;ii<n;ii++) //对A_temp赋值for(jj=0;jj<n;jj++)B(ii,jj)=X(jj,ii);for(ii=0;ii<n;ii++)L(ii,ii)=1.0;for (kk=0;kk<n;kk++) //总循环求出L U{ sum_answer=0;s_count=kk;for(pp=kk;pp<n;pp++)//doolitte 选主元{if (kk==0)sum_answer=0;elsefor(qq=0;qq<kk;qq++)sum_answer=sum_answer+L(pp,qq)*U(qq,pp);s=A_temp[pp][kk]-sum_answer;if (s>s_temp){s_temp=s;s_count=pp;}else if (s<(-1.0)*s_temp){s_temp=(-1.0)*s;s_count=pp;}else;}for(pp=0;pp<n;pp++) //换位A_temp{temp=A_temp[s_count][pp];A_temp[s_count][pp]=A_temp[kk][pp];A_temp[kk][pp]=temp;}for(pp=0;pp<kk;pp++) //换位L{temp=L(s_count,pp);L(s_count,pp)=L(kk,pp);L(kk,pp)=temp;}for(pp=0;pp<n;pp++) //换位B /////////////{temp=B(s_count,pp);B(s_count,pp)=B(kk,pp);B(kk,pp)=temp;}for (ii=kk;ii<n;ii++) //先求一行U{if(kk==0)U(kk,ii)=A_temp[kk][ii];else{sum_answer=0;for(qq=0;qq<kk;qq++)sum_answer=sum_answer+L(kk,qq)*U(qq,ii);U(kk,ii)=A_temp[kk][ii]-sum_answer;}}for(ii=kk+1;ii<n;ii++) //求出一列L{sum_answer=0;for(qq=0;qq<kk;qq++)sum_answer=sum_answer+L(ii,qq)*U(qq,kk);L(ii,kk)=(A_temp[ii][kk]-sum_answer)/U(kk,kk);}}//迭代Ly=bfor(kk=0;kk<n;kk++){Y[0][kk]=B(0,kk);for(ii=1;ii<n;ii++){sum_answer=0;for(jj=0;jj<ii;jj++){sum_answer=sum_answer+L(ii,jj)*Y[jj][kk];}Y[ii][kk]=B(ii,kk)-sum_answer;}//迭代Ux=yX((n-1),kk)=Y[n-1][kk]/U((n-1),(n-1));for(ii=n-2;ii>=0;ii--){sum_answer=0;for(jj=ii+1;jj<n;jj++){sum_answer=sum_answer+U(ii,jj)*X(jj,kk);}X(ii,kk)=(Y[ii][kk]-sum_answer)/U(ii,ii);}}for(ii=0;ii<n;ii++) // 给C赋值，直接转置for(jj=0;jj<n;jj++)C(ii,jj)=X(jj,ii);}。