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

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

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

3、报告版式要求打印版面要求：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 算法说明： -7 -1.2 ,算法程序： -8 -1.3 ,运行结果 -10 -2 .高斯列主元法解线性方程组 -10 -2.1 算法说明： -11 -2.2 算法程序： -11 -

5、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.2 算法程序： -25 -5.3 运行结果： -26 -6 .M次多项式曲线拟合 -27 -6.1 算法说明 -27 -6.2 算法程序： -28 -7 .不动点法解非线性方程 -33 -7.1 算法说明 -33 -7.3 运行结果：

6、-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*(

8、 K1+ 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

9、);long 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;/

10、h是步长 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(T

11、k+h,xk+h*f2,yk+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<<

12、setw(12)<<xk<<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;)陕西科技大学1.3,运行结果-10 -2.高斯列主元法解线性方程组2.1算法说明:首先

13、将线性方程组做成增广矩阵,对增广矩阵进行行变换.对于元素比, 在第i列中,第i行及以下的元素选取绝对值最大的元素,将该最大元素所 在的行与第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

14、"/#include "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(num

15、ber=0)包含名为enter的头函数 包含名为remove的头函数 包含名为judge的头函数调用enter函数 调用remove函数 调用judge函数陕西科技大学cout<<" 结果"<<endl;cout<<"此方程组无解"<<endl; /无解时予以提示)if(number=-1)(cout<<" 结果"<<endl;cout<<"此方程组有无穷多解"<<endl; / 有无穷多解时予以提示) if(numbe

16、r=1) ( cout<<" 结果"<<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&quo

17、t;<<i+1<<"="<<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 函数/寻找列主元,并移动

18、该行,最后化为上三角矩阵void remove(double *q,int m,int n)(int min=m;-12 -数值计算课程设计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;/通过换行以保证主对角线上的

19、元素是其所在位置及以下元素中绝对值最大的一个int t=0;for(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,in

20、t m,int n)(int r1=m;for(int i=0;i<m;i+)-13 -陕西科技大学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;)/求出增广矩

21、阵的秩求出系数矩阵的秩据不同的秩的情况返回不同的值2.3运行结果:-14 -数值计算课程设计-15 -陕西科技大学3.牛顿法解非线性方程组3.1算法说明设R.第1步：计算函数fi(Pk,qk)f2(Pk,qk)P重复上述过程.F(Pk)第2步：计算雅可比矩阵J(Pk)第3步：求线性方程组J(R) P的解P.第4步：计算下一点Ri Rfi( Pk,qk) fi(Pk,qk) xyf2(Pk,qk) f2(Pk,qk) xyF(R)3.2算法程序:#include<iostream>#include<cmath>#define N 3/非线性方程组中方程个数、未知量个数#d

22、efine EPsilon 0.0001 / 差向量1范数的上限#define Max 3/最大迭代次数using namesPace std;const int N2=2*N;int main()(void ff(float xxN,float yyN);/计算向量函数的因变量向量 yyN-16 -数值计算课程设计void ffjacobian(float xxN,float yyNN);/计算雅克比矩阵 yyNNvoid inv_jacobian(float yyNN,float invNN);/计算雅克比矩阵的逆矩阵 invvoid newdundiedai(float x0N, flo

23、at invNN,float y0N,float x1N);/由近似解向量x0计算近似解向量x1float x0N=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<<" "

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

25、;/计算差向量的1范数errornormerrornorm=0;for (i=0;i<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)-17 -陕西科技大学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*

26、y+z*z+z-6;cout<<" 向量函数的因变量向量是： "<<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

27、;yy12=2*z;yy20=2*x;yy21=2*y;yy22=2*z+1;cout<<"雅克比矩阵是："<<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;-18 -数值计算课程设计int i,j,k;cout<<"开始计算雅克比矩阵的逆矩

28、阵："<<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;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

29、;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=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<&l

30、t;""cout<<endl;cout<<endl;for (i=N-1;i>=0;i-)for(j=N2-1;j>=0;j-)augij=augij/augii;-19 -陕西科技大学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

31、=0;i<N;i+) for(j=0;j<N;j+)cout«invij«"cout«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

32、;i+)cout«x1i«"cout«endl;cout«endl;)3.3运行结果：-20-数值计算课程设计M /新M立堂关vt和逢试法De匕吧1牛迭代法用Ki初蛤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结果为：x=-1.36628;y=-0.366281;z=2.36628-21 -陕西科技大学4.龙贝格

33、求积分算法4.1 算法说明生成J K的逼近表R(J,K),并以R(J 1,J 1)为最终解来逼近积分bf (x)dx R(J, J)a逼近R(J,K)存在于一个特别的下三角矩阵中,第0列元素R(J,0)用基于2Jja,b子区间的连续梯形方法计算,然后利用龙贝格公式计算R(J,K).当1 K J时,第J行的元素为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&l

34、t;iostream>#include<cmath>#include<iomanip> using 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;cout<<"请输入两个端点a,b:"cin>>a>>

35、b;h=b-a;R00=h*(feval(a)+feval(b)while(err>0.0001)&&(J<N)|(J<4) /-22 -数值计算课程设计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<<"该矩阵为："<&

36、lt;endl;for(int i=0;i<=J;i+)for(int j=0;j<=J;j+)(cout<<setiosflags(ios:fixed)<<setprecision(8)<<setw(12)<<Rij; cout<<endl;cout<<"该式子在<<a<<","<<b<<"内的积分为"<<RJJ;cout<<endl;return 0; long double feval

37、(long double x) /函数的定义long double f;f=(x*x+x+1)*cos(x); return(f);-23 -陕西科技大学4.3运行结果:5.三次紧压样条插值5.1算法说明:-24 -数值计算课程设计三次紧压样条,确定S(%), S(Xn)(如果导数,这是“最正确选择),m0(d0S (x0)homN -一(S(Xn)5.2算法程序:#include<iostream>#include<iomanip>using namespace std;int main()const int n=4;const int N=n-1;/3long do

38、uble 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<<"求解经过点(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=

39、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=Hi4+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;

40、Bk=Bk-temp*Ck-1;Uk=Uk-temp*Uk-1;MN-1=UN-2/BN-2;-25 -陕西科技大学for(int k1=N-2;k1>0;k1-)Mk1=(Uk1-1-Ck1-1*Mk1+1)/Bk1-1;M0=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=:"<<e

41、ndl;for(int j=0;j<N;j+)for(int j1=0;j1<4;j1+)cout<<setw(12)<<Sjj1;cout<<endl;return 0;5.3运行结果:-26 -数值计算课程设计6,M次多项式曲线拟合6.1算法说明设xk» :有N个点,横坐标是确定的.最小二乘抛物线的系数表示为y f (x) Ax2 Bx C求解A, B和C的线性方程组为Nx： Ak 1Nx3 axk Ak1Nx2 Axk Ak1Nx； Bk1Nxk2 Bk1Nxk Bk1NN22xk Cykxkk1k1NNxk C ykxkk1k1

42、NNC ykk1-27 -陕西科技大学6.2算法程序:#include<iostream>#include<cmath>#include<iomanip>using namespace 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+

43、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);/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<<

44、;setw(13)<<F1p1k1;cout<<endl;cout<<endl;/输出Fcout<<"F=:"<<endl;for(int p2=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

45、;i2<m+1;i2+) for(int j2=0;j2<n;j2+)for(int k1=0;k1<m+1;k1+)F2j2k1=F1i2j2*Fj2k1;/每次输出F2cout<<"第"<<i2+1<<"次乘法："<<endl;for(int p3=0;p3<n;p3+)for(int k3=0;k3<m+1;k3+)-28 -数值计算课程设计cout<<setw(13)<<F2p3k3;cout<<endl; cout<<e

46、ndl;/得到矩阵Afor(int j3=0;j3<m+1;j3+) Ai2j3=0;/必须给此数赋初值for(int i3=0;i3<n;i3+) Ai2j3=Ai2j3+F2i3j3;/输出Acout<<"A=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>"<<end

47、l;for(int i4=0;i4<m+1;i4+) for(int j4=0;j4<n;j4+)Y1j4=F1i4j4*Yj4;/每次输出Y1cout<<"第"<<i4+1<<"次乘法："<<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<<

48、"B=F'*Y'=:"<<endl;for(int i6=0;i6<n;i6+)cout<<setw(13)<<Bi6;cout<<endl;invm+1m+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<<

49、"<以下计算矩阵C>"<<endl; for(int qq=0;qq<m+1;qq+) for(int w=0;w<n;w+)-29 -陕西科技大学B1w=invqqw*Bw;/每次输出B1cout<<"第"<<qq+1<<"次乘法："<<endl;for(int ww=0;ww<n;ww+)cout<<setw(13)<<B1ww;cout<<endl;/得到矩阵CCqq=0;/必须给此数赋初值for(int

50、 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;cout<<endl;return 0;long double inv_jacobian(long double yyNN,10ng double invNN)long double augNN2,L;int i,j,k;cout<<"开始计算A矩阵的逆矩阵："<<end

51、l;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<<"第零次消去："<<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=-au

52、gki/augii;for(j=i;j<N2;j+)augkj=augkj+L*augij;-30 -数值计算课程设计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 (k=i-1 ;k>=0;k-)L=-augki/augii;for(j=N2-1;j>=0;j-) augkj=augkj+L*augij； )cout

53、1;"第二次消去："«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(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);)-31 -陕西科技大学6.3运行结果:n=F1 *F=：G12.153S.B725151.215bSQ.2443245.2212.159S.S72E151.SIG3245-221E892.E38.0725151.215B8B.2443245.2215892.5787297151.215680.2443245.22158?2_573?2? . 739ZB98&80.2443245.2215892.57872?.73930081.9S632e+的63245.221E8T