典型数值算法的说明C++语言程序的设计开发

1、倏4猎菽人牵数值方法课程设计说明书典型数值算法的C+语言程序设计学生姓名：晏瑞学 号：200912010127院（系）：理学院专 业：数学与应用数学 091班指导教师：刘海峰2011 年6月15日陕西科技大学数值计算课程设计任务书理学院应用数学专业数学091班级学生：晏瑞 题目：典型数值算法的C+语言程序设计课程设计从 2011年5月20日起到 2011年6月25日1、课程设计的内容和要求（包括原始数据、技术要求、工作要求等）：每人需作10个算法的程序、必做6题、自选4题。对每个算法要求用C+语言进行编程。必选题：1、经典四阶龙格库塔法解一阶微分方程组2、高斯列主元法解线性方程组3、牛顿法解非

2、线性方程组4、龙贝格求积分算法5、三次样条插值算法（压紧样条）用 C+语言进行编程计算依据计算2S果，用Matlab画图并观察三次样条插值效果。6、M次多项式曲线拟合，据计算结果，用 Matlab画图并观察拟合效果。自选题：自选4道其他数值算法题目.每道题目重选次数不得超过5次.2、对课程设计成果的要求包括图表、实物等硬件要求：1）提交课程设计报告按照算法要求，用 C+语言设计和开发应用程序，提交由算法说明；程序 设计说明；系统技术文档（包括系统各模块主要流程图，软件测试方案与测试记 录、软件调试和修改记录、测试结论、运行情况记录），系统使用说明书，源程序代码为附录构成的课程设计报告。陕西科技

3、大学2）课程设计报告版式要求打印版面要求：A4纸，页边距：上 2cm,下2cm,左2.5cm、右2cm；字 体：正文宋体、小四号；行距：固定值 20;页眉1.5cm ,页脚1.75cm；页码位 于页脚居中打印；奇数页页眉“数值计算课程设计”，偶数页页眉“算法名称”, 页眉宋体小5号；段落及层次要求：每节标题以四号黑体左起打印（段前段后各 0.5行），节下为小节，以小四号黑体左起打印（段前段后各0.5行）。换行后以小四号宋体打印正文。节、小节分别以1、1.1、1.1.1依次标出，空一字符后接各部分的标题。当论文结构复杂，小节以下的标题，左起顶格书写，编号依次用（1）、（2） 或1）、2）顺序表示

4、。字体为小四号宋体。对条文内容采用分行并叙时，其编号用（a）、（b）或a）、b）顺序 表示，如果编号及其后内容新起一个段落，则编号前空两个中文字符。3）设计报告装订顺序与规范封面数值计算课程设计任务书目录数值计算设计课程设计报告正文设计体会及今后的改进意见参考文献（资料）左边缘装订3、课程设计工作进度计划：时间设计任务及要求第16周编写和调试程序并按要求撰写设计报告指导教师：日期：教研室主任：日期：目录1 .经典四阶龙格库塔法解一阶微分方程 -7 -1.1 算法说明： -71.2 ,算法程序： -81.3 ,运行结果 -10-2 .高斯列主元法解线性方程组 -10-2.1 算法说明： -11-

5、2.2 算法程序： -11 -2.3 运行结果： -14-3 .牛顿法解非线性方程组 -16-3.1 算法说明 -16-3.2 算法程序： -16-3.3 运行结果： -20-4 .龙贝格求积分算法 -22 -4.1 算法说明 -22-4.2 算法程序： -22-4.3 运彳丁结果： -24 -5 .三次紧压样条插值 -24 -5.1 算法说明： -24-5-数值计算课程设计5.2 算法程序： - 25 -5.3 运行结果： - 26 -6 .M 次多项式曲线拟合 - 27 - # -6.1 算法说明 -27-6.2 算法程序： -28-7 . 不动点法解非线性方程 -33-7.1 算法说明

6、-33-7.3 运行结果： -35-8 . 二分法解非线性方程 -36-8.1 算法说明： -36-8.4 运行结果： -38-9，龙格 -库塔法解微分方程 -39-9.1 算法说明： -39-9.2 算法程序： -40-9.3 运行界面： -41-10，递归梯形公式 -42-10.1 ，算法说明： -42-10.2 ，算法程序： -42-10.3 运行结果： -43-1.经典四阶龙格库塔法解一阶微分方程1.1 算法说明：龙格-库塔(Runge-Kutta) 方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基

7、础之上的。对于一精度的欧拉公式有：yi+1=yi+h*K1K1=f(xi,yi)当用点 xi 处的斜率近似值K1 与右端点xi+1 处的斜率K2 的算术平均值作为平均斜率K* 的近似值，那么就会得到二阶精度的改进欧拉公式：yi+1=yi+h*( K1+ K2)/2K1=f(xi,yi)K2=f(xi+h,yi+h*K1)依次类推，如果在区间xi,xi+1 内多预估几个点上的斜率值K1 、K2、Km,并用他们的加权平均数作为平均斜率K*的近似值，显然能构造出具有很高精度的高阶计算公式。经数学推导、求解， 可也就是在工程中应用广泛的经典龙格数值计算课程设计库塔算法：yi+1=yi+h*( K1+

8、2*K2 +2*K3+ K4)/6K1=f(xi,yi)K2=f(xi+h/2,yi+h*K1/2)K3=f(xi+h/2,yi+h*K2/2)K4=f(xi+h,yi+h*K3)通常所说的龙格-库塔法是指四阶而言的，我们可以仿二阶、三阶的情形推导出常用的标准四阶龙格-库塔法公式1.2 ,算法程序：#include<iostream>#include<iomanip>#define M 10using namespace std;int main()long double feval(long double ,long double ,long double );lon

9、g double feval1(long double ,long double ,long double );long double f4,g4;long double h,a,b,xa,ya;long double xM+1,yM+1,TM+1;cout<<"请输入区间左右端点a,b:"cin>>a>>b;cout<<endl;/xa,ya 是初值cout<<"请输入函数在左端点的初值xa,ya:"cin>>xa>>ya;cout<<endl;/h是步长

10、h=(b-a)/M;x0=xa;y0=ya;for(int i=0;i<=M;i+)Ti=a+h*i; / 给 T 赋值for(int k=0;k<=M;k+)f0=feval(Tk,xk,yk);g0=feval1(Tk,xk,yk);f1=feval(Tk+h/2,xk+h/2*f0,yk+h/2*g0);g1=feval1(Tk+h/2,xk+h/2*f0,yk+h/2*g0);f2=feval(Tk+h/2,xk+h/2*f1,yk+h/2*g1);g2=feval1(Tk+h/2,xk+h/2*f1,yk+h/2*g1);f3=feval(Tk+h,xk+h*f2,yk+

11、h*g2);g3=feval1(Tk+h,xk+h*f2,yk+h*g2);xk+1=xk+h/6*(f0+2*f1+2*f2+f3);yk+1=yk+h/6*(g0+2*g1+2*g2+f3);cout<<setw(12)<<"T"<<setw(12)<<"x"<<setw(12)<<"y"<<endl;for(k=0;k<=M;k+)cout<<setw(12)<<Tk<<setw(12)<<x

12、k<<setw(12)<<yk<<endl;return 0;long double feval(long double t,long double x,long double y)long double f;f=x+2*y;return f;long double feval1(long double t,long double x,long double y)long double f;f=3*x+2*y;return f;- 11 -1.3,运行结果2.1 算法说明：2.高斯列主元法解线性方程组首先将线性方程组做成增广矩阵，对增广矩阵进行行变换。对于元

13、素ai ,i ，在第 i 列中，第i 行及以下的元素选取绝对值最大的元素，将该最大元素所在的行与第i 行交换 ， 然后采用高斯消元法用新得到的ai ,i 所在的第i 行消去第 i 行以下的元素。依次进行直到an,n 。从而得到上三角矩阵。2.2 算法程序：本程序包含enter.h ， remove.h， judge.h及 main.cpp 这四个函数1， main.cpp 函数#include<iostream.h>#include<math.h>#include "enter.h"/#include "remove.h"/#in

14、clude "judge.h"/int main()int row,col;cout<<" 请输入方程的个数:"cin>>row;cout<<" 请输入未知数的个数加一:"cin>>col;double a100100,*p1,*p100;for(int i=0;i<row;i+)pi=ai;p1=p;enter(p1,row,col);/remove(p1,row,col);/int number=judge(p1,row,col); / if(number=0)包含名为ente

15、r 的头函数包含名为remove的头函数包含名为judge 的头函数调用enter 函数调用remove函数调用judge 函数数值计算课程设计cout<<" 结果 "<<endl;cout<<" 此方程组无解"<<endl; / 无解时予以提示if(number=-1)cout<<" 结果 "<<endl;cout<<" 此方程组有无穷多解"<<endl; / 有无穷多解时予以提示if(number=1)cout<

16、;<" 结果 "<<endl;cout<<" 此方程组有唯一解"<<endl;/有唯一解时，解出方程组的解double x100;xcol-2=arow-1col-1/arow-1row-1; /回代法求解上三角矩阵for(i=col-3;i>=0;i-) double t=aicol-1;for(int j=col-2;j>i;j-) t=t-aij*xj;xi=t/aii;for(i=0;i<col-1;i+)cout<<"x"<<i+1<&

17、lt;"="<<xi<<"t" / 依次输出方程组的解cout<<endl;return 0;2， enter.h 函数/ 输入系数矩阵void enter(double *q,int m,int n)cout<<" 请按行输入未知数前面的系数和等式右边的常数："<<endl;for(int i=0;i<m;i+)for(int j=0;j<n;j+) cin>>qij;3， remove.h 函数/寻找列主元，并移动该行，最后化为上三角矩阵void

18、remove(double *q,int m,int n)int min=m;if(n-1<m)min=n-1;for(int i=0;i<min-1;i+)int k=i;double max=qii;for(int j=i+1;j<m;j+)if(fabs(qji)>fabs(max) max=qji;k=j;/找到第i 列从 aii 开始的绝对值最大的元素if(k!=i)for(int j=0;j<n;j+)double mat=qij;qij=qkj;qkj=mat;/通过换行以保证主对角线上的元素是其所在位置及以下元素中绝对值最大的一个int t=0;f

19、or(j=i;j<m;j+)if(qji=0)t+;/在消元前判断aii 及其所在列以下元素是否都为零，不都为零再进行消元if(t!=m-i)for(int j=i+1;j<m;j+)for(int k=n-1;k>=i;k-) qjk=qjk-qik*qji/qii;for(i=0;i<n-2;i+)for(int j=i+1;j<m;j+)qji=0;/保证经消元后得到上三角阵4， judge.h 函数/用增广矩阵的秩判断该矩阵所对应的方程组的解的情况int judge(double *q,int m,int n)int r1=m;for(int i=0;i&

20、lt;m;i+)- 15 -数值计算课程设计-19 -int g=0;for(int j=0;j<n;j+)if(qij=0)g+;if(g=n)r1=r1-1;)int r2=m;for(i=0;i<m;i+)(int h=0;for(int j=0;j<n-1;j+)if(qij=0)h+;if(h=n-1)r2=r2-1;)if(r1=n-1&&r2=n-1)return 1;else if(r1=r2&&r1<n-1)return -1;else return 0;)/求出增广矩阵的秩求出系数矩阵的秩据不同的秩的情况返回不同的值2

21、.3运行结果:3.牛顿法解非线性方程组3.1算法说明设R已知。第1步：计算函数fi(Pk,qk)f2(Pk,qk)P重复上述过程。F(Pk)第2步：计算雅可比矩阵J(R)第3步：求线性方程组J(Pk) P的解P。第4步：计算下一点Ri Rfi(Pk,qk) fi(Pk,qk) xyf2(Pk,qk) f2(Pk,qk) xyF(Pk)3.2算法程序:#include<iostream>#include<cmath>#define N 3/非线性方程组中方程个数、未知量个数#define Epsilon 0.0001 / 差向量1范数的上限#define Max 3/最大

22、迭代次数using namespace std;const int N2=2*N;int main()void ff(float xxN,float yyN);/计算向量函数的因变量向量 yyN陕西科技大学void ffjacobian(float xxN,float yyNN);/计算雅克比矩阵yyNNvoid inv_jacobian(float yyNN,float invNN);/计算雅克比矩阵的逆矩阵invvoid newdundiedai(float x0N, float invNN,float y0N,float x1N);/由近似解向量x0 计算近似解向量x1float x0N

23、=0,0,0,y0N,jacobianNN,invjacobianNN,x1N,errornorm;int i,j,iter=0;/ 如果取消对x0 的初始化，撤销下面两行的注释符，就可以由键盘向x0 读入初始近似解向量/for( i=0;i<N;i+)/cin>>x0i;cout<<" 初始近似解向量："<<endl;for (i=0;i<N;i+)cout<<x0i<<" " cout<<endl;cout<<endl;doiter=iter+1;cout

24、<<" 第 "<<iter<<" 次迭代开始"<<endl;/ 计算向量函数的因变量向量y0ff(x0,y0);/ 计算雅克比矩阵jacobianffjacobian(x0,jacobian);/ 计算雅克比矩阵的逆矩阵invjacobianinv_jacobian(jacobian,invjacobian);/ 由近似解向量x0 计算近似解向量x1newdundiedai(x0, invjacobian,y0,x1);/ 计算差向量的1 范数 errornormerrornorm=0;for (i=0;i

25、<N;i+)errornorm=errornorm+fabs(x1i-x0i);if (errornorm<Epsilon) break;for (i=0;i<N;i+) x0i=x1i; while (iter<Max);return 0;void ff(float xxN,float yyN)float x,y,z;int i;x=xx0;y=xx1;z=xx2;yy0=x*x-x+y*y+z*z-5;yy1=x*x+y*y-y+z*z-4;yy2=x*x+y*y+z*z+z-6;cout<<" 向量函数的因变量向量是："<&l

26、t;endl;for( i=0;i<N;i+)cout<<yyi<<" "cout<<endl;cout<<endl;void ffjacobian(float xxN,float yyNN)float x,y,z;int i,j;x=xx0;y=xx1;z=xx2;/jacobian have n*n elementyy00=2*x-1;yy01=2*y;yy02=2*z;yy10=2*x;yy11=2*y-1;yy12=2*z;yy20=2*x;yy21=2*y;yy22=2*z+1;cout<<&quo

27、t; 雅克比矩阵是："<<endl;for( i=0;i<N;i+)for(j=0;j<N;j+) cout<<yyij<<" "cout<<endl;cout<<endl;void inv_jacobian(float yyNN,float invNN)float augNN2,L;- # -陕西科技大学int i,j,k;cout<<"开始计算雅克比矩阵的逆矩阵："<<endl;for (i=0;i<N;i+) for(j=0;j<N

28、;j+)augij=yyij;for(j=N;j<N2;j+)if(j=i+N) augij=1;else augij=0;for (i=0;i<N;i+) for(j=0;j<N2;j+)cout<<augij<<""cout<<endl;cout<<endl;for (i=0;i<N;i+)for (k=i+1;k<N;k+)L=-augki/augii;for(j=i;j<N2;j+) augkj=augkj+L*augij; for (i=0;i<N;i+) for(j=0;j

29、<N2;j+)cout<<augij<<""cout<<endl;cout<<endl;for (i=N-1;i>0;i-)for (k=i-1;k>=0;k-)L=-augki/augii;for(j=N2-1;j>=0;j-) augkj=augkj+L*augij; for (i=0;i<N;i+) for(j=0;j<N2;j+)cout<<augij<<""cout<<endl;cout<<endl;for (i=

30、N-1;i>=0;i-)for(j=N2-1;j>=0;j-)augij=augij/augii;for (i=0;i<N;i+) for(j=0;j<N2;j+)cout«augij«"cout«endl;for(j=N;j<N2;j+)invij-N=augij；)cout«endl;cout«"雅克比矩阵的逆矩阵："«endl;for (i=0;i<N;i+) for(j=0;j<N;j+)cout«invij«"cout

31、71;endl;)cout«endl;)void newdundiedai(float xON, float invNN,float yON,float x1N) (int ij;float sum=0;for(i=0;i<N;i+)sum=0;for(j=0;j<N;j+)sum=sum+invij*yOj;x1i=x0i-sum;)cout«"近似解向量："«endl;for (i=0;i<N;i+)cout«x1i«"cout«endl;cout«endl;)3.3运行结

32、果：-#-数值计算课程设计M /新M立堂关vt和逢试法De匕吧1牛迭代法用Ki结果为：x=-1.36628;y=-0.366281;z=2.36628初蛤Sffl晌量tJ M 0黯皤感%向量是:5T养克比矩阵是； -10 6 0 -1 G 0 0 1卜力铝+算雅克比矩阵的逆矩阵: -1 QQ10a0 -1s01a0 Q 1 0 0 1-i Bsi0a0 -1q01ae n i s 0 1-1 0fi100(J Yg(d10-25 -数值计算课程设计4.龙贝格求积分算法4.1 算法说明生成J K的逼近表R(J,K),并以R(J 1,J 1)为最终解来逼近积分bf (x)dx R(J, J)a逼近

33、R(J,K)存在于一个特别的下三角矩阵中，第0列元素R(J,0)用基于2J个a,b子区间的连续梯形方法计算，然后利用龙贝格公式计算R(J,K)。当1 K J时，第J行的元素为-31 -R(J,K) R(J,K 1)R(J,K 1) R(J 1,K 1)声明在指定函数的指定点处的函数函数的循环当|R(J,J) R(J 1,J 1)| tol时，程序在第(J 1)行结束。4.2 算法程序:/龙贝格方法求定积分/函数的修改可以在feval函数的定义中修改#include<iostream>#include<cmath>#include<iomanip> using

34、 namespace std;#define N 4 /定义全局变量 Nint main()long double feval(long double );/long double a,b,h;long double R1010=0;long double x;int M=1;long double err=1;int J=0;long double s=0;coutvv"请输入两个端点a,b:"cin>>a>>b;h=b-a;R00=h*(feval(a)+feval(b)while(err>0.0001)&&(J<N)

35、|(J<4) /J=J+1; h=h s=0;for(int p=1;p<=M;p+) (x=a+h*(2*p-1);s=s+feval(x);RJ0=RJ-10/2+h*s;M=2*M;for(int K=1;K<=J;K+) (RJK=RJK-1+(RJK-1-RJ-1K-1)/(4AK-1); err=fabs(RJ-1J-1-RJK); cout<<"该矩阵为："<<endl;for(int i=0;i<=J;i+)for(int j=0;j<=J;j+)(cout<<setiosflags(ios:

36、fixed)<<setprecision(8)<<setw(12)<<Rij; cout<<endl;cout<<"该式子在<<a<<","<<b<<"内的积分为"<<RJJ;cout<<endl;return 0; long double feval(long double x) /函数的定义long double f;f=(x*x+x+1)*cos(x); return(f);4.3运行结果:5.三次紧压样条插

37、值5.1算法说明:三次紧压样条，确定SU), S(%)(如果导数已知，这是“最佳选择”),m0(d0 S (x0)h0mN -一(S(Xn)5.2算法程序:#include<iostream>#include<iomanip>using namespace std;int main()const int n=4;const int N=n-1;/3long double Xn=0,1,2,3,Yn=0,0.5,2.0,1.5,dx0=0.2,dxn=-1,Hn-1,H1n-1,Dn-1,Un-2,AN-2,BN-1,CN-1,MN+1,temp,SN4;cout<

38、<"求解经过点(0, 0.0),(1 , 0.5),(2 , 2.0)和(3, 1.5),而且一阶导数边界条件S'(0)=0.2和S'(3)=-1的三次压紧样条曲线："；for(int i=0;i<n-1;i+)Hi=Xi+1-Xi;for(int i1=0;i1<n-1;i1+)H1i1=Yi1+1-Yi1;for(int i2=0;i2<n-1;i2+)Di2=H1i2/Hi2;for(int i3=0;i3<n-1;i3+)Ui3=6*(Di3+1-Di3);for(int i4=0;i4<N-2;i4+)Ai4=H

39、i4+1;for(int i5=0;i5<N-1;i5+)Bi5=2*(Hi5+Hi5+1);for(int i6=0;i6<N-1;i6+)Ci6=Hi6+1;B0=B0-H0/2;U0=U0-3*(D0-dx0);BN-2=BN-2-HN-1/2;UN-2=UN-2-3*(dxn-DN-1);for(int k=1;k<N-1;k+)temp=Ak-1/Bk-1;Bk=Bk-temp*Ck-1;Uk=Uk-temp*Uk-1;MN-1=UN-2/BN-2;for(int k1=N-2;k1>0;k1-)Mk1=(Uk1-1-Ck1-1*Mk1+1)/Bk1-1;M0

40、=3*(D0-dx0)/H0-M1/2;MN=3*(dxn-DN-1)/HN-1-MN-1/2;for(int k2=0;k2<N;k2+)Sk20=(Mk2+1-Mk2)/(6*Hk2);Sk21=Mk2/2;Sk22=Dk2-Hk2*(2*Mk2+Mk2+1)/6;Sk23=Yk2;cout<<"S=:"<<endl;for(int j=0;j<N;j+)for(int j1=0;j1<4;j1+)cout<<setw(12)<<S皿1;cout<<endl;return 0;5.3运行结果:

41、数值计算课程设计6.M 次多项式曲线拟合6.1 算法说明设（xk,yk） N1有N个点，横坐标是确定的。最小二乘抛物线的系数表示为y f (x) Ax2 Bx C求解A, B和C的线性方程组为- 33 -Nxk4 A k1Nx3 Axk A k1Nx2 Axk A k1Nxk3 B k1Nxk2 B k1Nxk B k1NN22xk C ykxkk1k1NNxk Cykxkk1k1NNC ykk1数值计算课程设计6.2算法程序:#include<iostream>#include<cmath>#include<iomanip>using namespace

42、std;const int n=6,m=5,N=6,N2=2*N;int main()long double inv_jacobian(long double yyNN,10ng double invNN); long double Xn=0.25,1.0,1.5,2.0,2.4,5.0,Yn=23.1,1.68,1.0,0.84,0.826,1.2576,Bm+1,Fnm+1=0,F1m+1n,F2nm+1,Am+1m+1,Y1n,invm+1m+1,Cm+1,B1m+1;for(int i=0;i<m+1;i+)for(int j=0;j<n;j+)Fji=pow(Xj,i);

43、/F的转置for(int i1=0;i1<m+1;i1+)for(int j1=0;j1<n;j1+)F1i1j1=Fj1i1;/输出F1cout<<"F 的转置 F'=:"<<endl;for(int p1=0;p1<m+1;p1+)for(int k1=0;k1<n;k1+)cout<<setw(13)<<F1p1k1;cout<<endl;cout<<endl;/输出Fcout<<"F=:"<<endl;for(int p

44、2=0;p2<n;p2+)for(int k2=0;k2<m+1;k2+)cout<<setw(13)<<Fp2k2;cout<<endl;cout<<endl;cout<<"<以下计算矩阵A>"<<endl;/F的转置F1乘以F得A矩阵for(int i2=0;i2<m+1;i2+) for(int j2=0;j2<n;j2+)for(int k1=0;k1<m+1;k1+)F2j2k1=F1i2j2*Fj2k1;/每次输出F2cout<<&quo

45、t;第"<<i2+1<<"次乘法："<<endl;for(int p3=0;p3<n;p3+)for(int k3=0;k3<m+1;k3+)cout<<setw(13)<<F2p3k3;cout<<endl; cout<<endl;/得到矩阵Afor(int j3=0;j3<m+1;j3+) Ai2j3=0;/必须给此数赋初值for(int i3=0;i3<n;i3+) Ai2j3=Ai2j3+F2i3j3;/输出Acout<<"A=

46、F'*F=:"<<endl;for(int p=0;p<m+1;p+)for(int k=0;k<m+1;k+)cout<<setw(13)<<Apk;cout<<endl;cout<<"<以下计算矩阵B>"<<endl;for(int i4=0;i4<m+1;i4+) for(int j4=0;j4<n;j4+)Y1j4=F1i4j4*Yj4;/每次输出Y1cout<<"第"<<i4+1<<&

47、quot;次乘法："<<endl;for(int p4=0;p4<n;p4+)cout<<setw(13)<<Y1p4;cout<<endl;/得到矩阵BBi4=0;/必须给此数赋初值for(int i5=0;i5<m+1;i5+)Bi4=Bi4+Y1i5;cout<<"B=F'*Y'=:"<<endl;for(int i6=0;i6<n;i6+)cout<<setw(13)<<Bi6;cout<<endl;invm+1m+

48、1=inv_jacobian(A,inv);cout<<"A矩阵的逆矩阵："<<endl;for (i=0;i<N;i+) for(int j0=0;j0<N;j0+)cout<<setw(13)<<invij0;cout<<endl;cout<<"<以下计算矩阵C>"<<endl; for(int qq=0;qq<m+1;qq+) for(int w=0;w<n;w+)B1w=invqqw*Bw;/每次输出B1cout<<

49、"第"<<qq+1<<"次乘法："<<endl;for(int ww=0;ww<n;ww+)cout<<setw(13)<<B1ww;cout<<endl;/得到矩阵CCqq=0;/必须给此数赋初值for(int e=0;e<m+1;e+)Cqq=Cqq+B1e;cout<<"C=A'*B="<<endl;for(int tt=0;tt<m+1;tt+)cout<<setw(13)<<Ctt

50、;cout<<endl;return 0;long double inv_jacobian(long double yyNN,10ng double invNN)long double augNN2,L;int i,j,k;cout<<"开始计算A矩阵的逆矩阵："<<endl;for (i=0;i<N;i+) for(j=0;j<N;j+)augij=yyij;for(j=N;j<N2;j+)if(j=i+N) augij=1;else augij=0;cout<<"第零次消去："<

51、<endl;for (i=0;i<N;i+) for(j=0;j<N2;j+)cout<<augij<<"/"cout<<endl;cout<<endl;for (i=0;i<N;i+)for (k=i+1;k<N;k+)L=-augki/augii;for(j=i;j<N2;j+)augkj=augkj+L*augij;-39 -cout«"第一次消去："«endl;for (i=0;i<N;i+) for(j=0;j<N2;j+) co

52、utvvaugijvvr cout«endl;)cout«endl;for (i=N-1;i>0;i-)(for (k=i-1 ;k>=0;k-)L=-augki/augii;for(j=N2-1;j>=0;j-) augkj=augkj+L*augij； )cout«"第二次消去："«endl;for (i=0;i<N;i+) for(j=0;j<N2;j+) coutvvaugijvvr cout«endl;)cout«endl;for (i=N-1;i>=0;i-) for

53、(j=N2-1;j>=0;j-)augiO=augij/augii;cout«"第三次消去："«endl;for (i=0;i<N;i+) for(j=0;j<N2;j+) coutvvaugijvvr cout«endl;for(j=N;j<N2;j+) invij-N=augij；) cout«endl;return(invm+1m+1);)6.3运行结果:n=F1 *F=：G12.153S.B725151.215bSQ.2443245.2212.159S.S72E151.SIG3245-221E892.E

54、38.0725151.215B8B.2443245.2215892.5787297151.215680.2443245.22158?2_573?2? . 739ZB98&80.2443245.2215892.57872?.73930081.9S632e+的63245.221E8T2.578729,7393Q81.95632e*00t?.?7385e*S06pl®慑矩酗23 .11.6g10.6261.2576酬欣乘法：5.7/1.661_51.60l.?8246,260情啾乘法：1-443751.6U2.25J.364-VS7763L.44用欷乘法： 0.36U93S1.6S3.3