线性方程组的直接解法 实验报告
解线性方程组的直接方法实验报告范本
解线性方程组的直接方法实验报告Record the situation and lessons learned, find out the existing problems andform future countermeasures.姓名：___________________单位：___________________时间：___________________编号：FS-DY-20246 解线性方程组的直接方法实验报告解线性方程组的直接方法实验报告1.实验目的：1、通过该课题的实验，体会模块化结构程序设计方法的优点；2、运用所学的计算方法，解决各类线性方程组的直接算法；3、提高分析和解决问题的能力，做到学以致用；4、通过三对角形线性方程组的解法，体会稀疏线性方程组解法的特点。
2.实验过程：实验代码：#include "stdio.h"#include "math.h"#includeusing namespace std;//Gauss法void lzy(double **a,double *b,int n) {int i,j,k;double l,x[10],temp;for(k=0;k=0;i--){temp=0;for(j=i+1;j=0;k--){temp=0;for(m=k+1;m0)d[i-1]=a[i][i-1];}a1[0]=a0[0];for(i=0;i=0;i--)x[i]=y[i]-b1[i]*x[i+1];for(i=0;i>n;cout<<" 请输入系数矩阵: ";for(i=0;i<n;i++)for(j=0;j<n;j++){Foonshion图文设计有限公司Fonshion Design Co., Ltd。
实验五 解线性方程组的直接法
毕节学院实验报告实验名称：解线性方程组的直接法实验报告序号： 5 组别姓名罗晟同组实验者实验项目解线性方程组的直接法实验日期2012年11月17日实验类别□√1、验证性实验或基础性实验；□2、综合性实验□3、设计性实验；□4、创新性实验和研究性实验；教师评语实验成绩指导教师（签名）赖志柱年月日实验目的掌握线性方程组直接解法的基本思想，进一步熟练掌握高斯消去法，提高编程能力和解算线性方程组问题的实践技能。
实验任务与要求用一种高级语言（推荐用MATLAB）编写高斯消去法具体实现的函数，要求输入参数包含方程组的系数矩阵及常数项向量，输出方程组的解，并使用若干个有代表性的例子进行调试。
小组合作分工说明：实验过程及内容：function [x,det,index] = LGuass(A,b)[n,m] = size(A);nb = length(b);if n ~= merror('矩阵A的行列必须相等！');return;endif m ~= nberror('b矩阵的个数必须与矩阵A的行数相等！');return;end%开始计算，先赋初值index = 1;det = 1;x = zeros(n,1);for k = 1:n-1%选主元a_max = 0;for i = k:nif abs(A(i,k)) > a_maxa_max = abs(A(i,k));r = i;endendif a_max < 1e-20index = 0;return;end%交换两行if r > kfor j = k:nz = A(k,j);A(k,j) = A(r,j);A(r,j) = z;endz = b(k);b(k) = b(r);b(r) = z;det = -det;end%消元过程for i = k+1:nm = A(i,k)/A(k,k);for j = k+1:nA(i,j) = A(i,j) - m*A(k,j);endb(i) = b(i) - m*b(k);enddet = det*A(k,k);enddet = det*A(n,n);%回代过程if abs(A(n,n)) < 1e-10index = 0;return;endfor k = n:-1:1for j = k+1:nb(k) = b(k) - A(k,j)*x(j);endx(k)=b(k)/A(k,k);end。
实验2 线性方程组的直接解法
实验2 线性方程组解的直接解法【实验目的】1、验证高斯消去法和三角分解法2、掌握直接求解线性方程组的常用算法：列主元高斯消去法、LU 分解法等3、记录运行结果，回答问题，完成实验报告【实验要求】根据题目要求，用Matlab 完成下列实验内容。
书写实验报告。
【实验内容】1、利用Matlab 求解线性方程组考虑下面给出的线性方程组形式为Ax =B 。
其中，A 和B 均为给定矩阵1111n1,n nn n a a b A B a a b ⎛⎫⎛⎫ ⎪ ⎪== ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭ 可以由给定的A 和B 矩阵构造出解的判定矩阵C （可以用C=[A B]命令来生成C ）11111n n nn n a a b C a a b ⎛⎫ ⎪= ⎪ ⎪⎝⎭ 这里先可以给出线性方程组有解的判定定理：（2）当rank(A )=rank(C )=r ＜n 时，方程组有无穷多解。
具体求解方法留给同学们实验课后自学。
（3）当rank(A )＜rank(C )时，方程组为矛盾方程，只能解出方程的最小二乘解。
这部分内容将在第四章学习。
【问题1】：求解下面给出的线性方程组。
记录所用命令和求解结果。
12345432141324341322x ⎛⎫⎛⎫ ⎪ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭ 【思考1】：Matlab 对线性方程组求解是否要考虑主元的选取？为什么？可以用下列例子试试 12120.000112x x x x +=⎧⎨+=⎩2、用Matlab 编写高斯消去法具体实现的函数要求输入参数包含方程组的系数矩阵及常数项向量，输出方程组的解，并使用问题1的例子进行调试。
高斯顺序消去法:设线性方程组AX ＝b 。
1.消元过程设0)(≠k kk a ,对1,,2,1-=n k 计算()()(1)()()(1)()()/,1,2,,k k ik ik kk k k k ij ij ik kj k k k i i ik k m a a a a m a i j k k n b b m b ++⎧=⎪=-=++⎨⎪=-⎩K 其中⒉回代过程()()()()()1/1,,2,1()/n n n n nn n i i i ii ij j ii j i x b a i n x b a x a =+⎧=⎪=-⎨=-⎪⎩∑K 其中。
线性方程组的直接解法实验报告
本科实验报告
课程名称：数值计算方法B
实验项目：线性方程组的直接解法
最小二乘拟合多项式
实验地点：ZSA401
专业班级：学号：201000
学生姓名：
指导教师：李志
2012年4月13日
for(i=1;i<=n;i++)
{
for(j=1;j<=n+1;j++)
printf("%lf\t",A[i][j]);
printf("\n");
}
double answer[N];
Gauss_eliminate(n,answer);
/*输出解*/
for(i=1;i<=n;i++)
printf("a[%d]=%lf\t",i-1,answer[i]);
getchar();
getchar();
}
四、实验结果与讨论、心得
讨论、心得:
刚开始调试代码的时候有时候就是很小的错误导致整个程序不能运行，需要我们一步一步慢慢来，经过无数次的检查程序错误的原因，以及在老师的帮助下，完成了这次实验。
这段时间的实验课提高了我的分析问题，解决问题的能力，特别提高了对一个程序的整。
解线性方程组的直接方法实验报告
解线性方程组的直接方法实验报告解线性方程组的直接方法实验报告1.实验目的：1、通过该课题的实验，体会模块化结构程序设计方法的优点；2、运用所学的计算方法，解决各类线性方程组的直接算法；3、提高分析和解决问题的能力，做到学以致用；4、通过三对角形线性方程组的解法，体会稀疏线性方程组解法的特点。
2.实验过程：实验代码：#include "stdio.h"#include "math.h"#includeusing namespace std;//Gauss法void lzy(double **a,double *b,int n){int i,j,k;double l,x[10],temp;for(k=0;k<n-1;k++){for(j=k,i=k;j<n;j++){if(j==k)temp=fabs(a[j][k]);else if(temp<fabs(a[j][k])) {temp=fabs(a[j][k]);i=j;}}if(temp==0){cout<<"无解" ; return;}else{for(j=k;j<n;j++){temp=a[k][j];a[k][j]=a[i][j];a[i][j]=temp;}temp=b[k];b[k]=b[i];b[i]=temp;}for(i=k+1;i<n;i++){l=a[i][k]/a[k][k];for(j=k;j<n;j++)a[i][j]=a[i][j]-l*a[k][j];b[i]=b[i]-l*b[k];}}if(a[n-1][n-1]==0){cout<<"无解" ; return;}x[n-1]=b[n-1]/a[n-1][n-1]; for(i=n-2;i>=0;i--){temp=0;for(j=i+1;j<n;j++)temp=temp+a[i][j]*x[j];x[i]=(b[i]-temp)/a[i][i];}for(i=0;i<n;i++){printf("x%d=%lf ",i+1,x[i]);printf(" ");}}//平方根法void pfg(double **a,double *b,int n) {int i,k,m;double x[8],y[8],temp;for(k=0;k<n;k++){temp=0;for(m=0;m<k;m++)temp=temp+pow(a[k][m],2);if(a[k][k]<temp)return;a[k][k]=pow((a[k][k]-temp),1.0/2.0);for(i=k+1;i<n;i++){temp=0;for(m=0;m<k;m++)temp=temp+a[i][m]*a[k][m]; a[i][k]=(a[i][k]-temp)/a[k][k]; }temp=0;for(m=0;m<k;m++)temp=temp+a[k][m]*y[m];y[k]=(b[k]-temp)/a[k][k];}x[n-1]=y[n-1]/a[n-1][n-1];for(k=n-2;k>=0;k--){temp=0;for(m=k+1;m<n;m++)temp=temp+a[m][k]*x[m];x[k]=(y[k]-temp)/a[k][k];}for(i=0;i<n;i++){printf("x%d=%lf ",i+1,x[i]);printf(" ");}}//追赶法void zgf(double **a,double *b,int n){int i;double a0[10],c[10],d[10],a1[10],b1[10],x[10],y[10]; for(i=0;i<n;i++) {a0[i]=a[i][i];if(i<n-1)c[i]=a[i][i+1];if(i>0)d[i-1]=a[i][i-1];}a1[0]=a0[0];for(i=0;i<n-1;i++){b1[i]=c[i]/a1[i];a1[i+1]=a0[i+1]-d[i+1]*b1[i];}y[0]=b[0]/a1[0];for(i=1;i<n;i++)y[i]=(b[i]-d[i]*y[i-1])/a1[i];x[n-1]=y[n-1];for(i=n-2;i>=0;i--)x[i]=y[i]-b1[i]*x[i+1];for(i=0;i<n;i++){printf("x%d=%lf ",i+1,x[i]);printf(" ");}}int main{int n,i,j;double **A,**B,**C,*B1,*B2,*B3;A=(double **)malloc(n*sizeof(double)); B=(double **)malloc(n*sizeof(double));C=(double **)malloc(n*sizeof(double));B1=(double *)malloc(n*sizeof(double));B2=(double *)malloc(n*sizeof(double));B3=(double *)malloc(n*sizeof(double));for(i=0;i<n;i++){A[i]=(double *)malloc((n)*sizeof(double)); B[i]=(double *)malloc((n)*sizeof(double)); C[i]=(double *)malloc((n)*sizeof(double)); } cout<<"第一题(Gauss列主元消去法):"<<endl<<endl; cout<<"请输入阶数n:"<<endl;cin>>n;cout<<" 请输入系数矩阵: ";for(i=0;i<n;i++)for(j=0;j<n;j++){。
实验五 解线性方程组的直接方法
实验五 解线性方程组的直接方法实验5.1 （主元的选取与算法的稳定性）问题提出：Gauss 消去法是我们在线性代数中已经熟悉的。
但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss 消去法作为数值算法的稳定性呢？Gauss 消去法从理论算法到数值算法，其关键是主元的选择。
主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。
实验内容：考虑线性方程组n n n R b R A b Ax ∈∈=⨯,,编制一个能自动选取主元，又能手动选取主元的求解线性方程组的Gauss 消去过程。
实验要求：（1）取矩阵⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=1415157,6816816816 b A ，则方程有解T x )1,,1,1(* =。
取n=10计算矩阵的条件数。
让程序自动选取主元(顺序消元)，结果如何？（2）现选择程序中手动选取主元的功能。
每步消去过程总选取按模最小或按模尽可能小的元素作为主元，观察并记录计算结果。
若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。
（3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。
（4）选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。
重复上述实验，观察记录并分析实验结果。
实验5.2（线性代数方程组的性态与条件数的估计）问题提出：理论上，线性代数方程组b Ax =的摄动满足⎪⎪⎭⎫ ⎝⎛∆+∆∆-≤∆-b b A A A A A cond x x 11)(矩阵的条件数确实是对矩阵病态性的刻画，但在实际应用中直接计算它显然不现实，因为计算1-A 通常要比求解方程b Ax =还困难。
实验内容：MATLAB 中提供有函数“condest ”可以用来估计矩阵的条件数，它给出的是按1-范数的条件数。
计算方法-解线性方程组的直接法实验报告
cout<<endl;
for(k=i+1;k<m;k++)
{
l[k][i]=a[k][i]/a[i][i];
for(r=i;r<m+1;r++) /*化成三角阵*/
a[k][r]=a[k][r]-l[k][i]*a[i][r];
}
}
x[m-1]=a[m-1][m]/a[m-1][m-1];
{
int i,j;
float t,s1,s2;
float y[100];
for(i=1;i<=n;i++) /*第一次回代过程开始*/
{
s1=0;
for(j=1;j<i;j++)
{
t=-l[i][j];
s1=s1+t*y[j];
}
y[i]=(b[i]+s1)/l[i][i];
}
for(i=n;i>=1;i--) /*第二次回代过程开始*/
s2=s2+l[i][k]*u[k][r];
l[i][r]=(a[i][r]-s2)/u[r][r];
}
}
printf("array L:\n");/*输出矩阵L*/ for(i=1;i<=n;i++)
{
for(j=1;j<=n;j++)
printf("%7.3f ",l[i][j]);
printf("\n");
{
s2=0;
for(j=n;j>i;j--)
求解线性方程组的直接方法算法实验报告
求解线性方程组的直接方法（2学时）一 实验目的1．掌握求解线性方程组的高斯消元法及列主元素法；2. 掌握求解线性方程组的克劳特法；3. 掌握求解线性方程组的平方根法。
二 实验内容1．用高斯消元法求解方程组（精度要求为610-=ε）：1231231233272212240x x x x x x x x x -+=⎧⎪-+-=-⎨⎪-+=⎩2．用克劳特法求解上述方程组（精度要求为610-=ε）。
3. 用平方根法求解上述方程组（精度要求为610-=ε）。
4. 用列主元素法求解方程组（精度要求为610-=ε）：1231231233432222325x x x x x x x x x -+=⎧⎪-+-=⎨⎪--=-⎩三 实验步骤（算法）与结果1. #include<stdio.h>main(){ float a11,a12,a13,b1,a21,a22,a23,b2,a31,a32,a33,b3,l21,l31,l32,u11,u12,u13,u22,u23,u33,z1,z2,z3,x1,x2,x3;printf("enter a11,a12,a13,b1,a21,a22,a23,b2,a31,a32,a33,b3 :"); scanf("%f%f%f%f%f%f%f%f%f%f%f%f/n",&a11,&a12,&a13,&b1,&a21,&a22, &a23,&b2,&a31,&a32,&a33,&b3 );l21=a12/a11;l31=a31/a11;u22=a22-l21*a12;l32=(a32-l31*a12)/u22;u23=a23-l21*a13;u33=a33-l31*a13-l32*u23;z2=b2-l21*b1;z3=b3-l31*b1-l32*z2;printf("u22=%fu23=%fz2=%fu33=%fz3=%f`````",u22,u23,z2,u33,z3); x3=z3/u33;x2=(z2-u23*x3)/u22;x1=(b1-a13*x3-a12*x2)/a11;printf("x1=%f x2=%f x3=%f ",x1,x2,x3);return 0;}2．#include<stdio.h>main(){float a11,a12,a13,b1,a21,a22,a23,b2,a31,a32,a33,b3,l22,l32,l33,u11,u12,u13,u22,u23,u33,z1,z2,z3,x1,x2,x3;printf("enter a11,a12,a13,b1,a21,a22,a23,b2,a31,a32,a33,b3 :");scanf("%f%f%f%f%f%f%f%f%f%f%f%f/n",&a11,&a12,&a13,&b1,&a21,&a22,&a23,&b2,&a31,&a32,&a33,&b3 );u11=1;u22=1; u33=1;u12=a12/a11;u13=a13/a11;z1=b1/a11;l22=a22-a21*u12;u23=(a23-a21*u13)/l22;z2=(b2-a21*z1)/l22;l32=a32-a31*u12;l33=a33-a31*u13-l32*u23;z3=(b3-a31*z1-l32*z2)/l33;printf("u11=%f u12=%f u13=%f z1=%f u22=%f u23=%f z2=%f u33=%f z3=%f------",u11,u12,u13,z1,u22,u23,z2,u33,z3);x3=z3;x2=z2-u23*x3;x1=z1-u13*x3-u12*x2;printf("x1=%f x2=%f x3=%f ",x1,x2,x3);getch();return 0;}3． #include<stdio.h>#include<math.h>{float a11,a12,a13,b1,a21,a22,a23,b2,a31,a32,a33,b3,l11,l12,l13,l23,l21,l22,l31,l32,l33,z1,z2,z3,x1,x2,x3;printf("enter a11,a12,a13,b1,a21,a22,a23,b2,a31,a32,a33,b3 :");scanf("%f%f%f%f%f%f%f%f%f%f%f%f/n",&a11,&a12,&a13,&b1,&a21,&a22, &a23,&b2,&a31,&a32,&a33,&b3 );l11=sqrt(a11);l21=a21/l11; l31=a31/l11;l22=sqrt(a22-l21*l21);l32=(a32-l21*l31)/l22;l33=sqrt(a33-l31*l31-l32*l32);z1=b1/l11;z2=(b2-l21*z1)/l22;z3=(b3-l31*z1-l32*z2)/l33;printf("l11=%f z1=%f l22=%f z2=%f l33=%fz3=%f---",l11,z1,l22,z2,l33,z3);x3=z3/l33;x2=(z2-l32*x3)/l22;x1=(z1-l31*x3-l21*x2)/l11;printf("x1=%f x2=%f x3=%f ",x1,x2,x3);getch();return 0;}4． #include "stdio.h"#include "math.h"main(){ float a11,a12,a13,b1,a21,a22,a23,b2,a31,a32,a33,b3,l21,l31,A22,A23,d1,A32,A33,d2,l32,a,d3,x1,x2,x3,A,B,C,D;printf("enter a11,a12,a13,b1,a21,a22,a23,b2,a31,a32,a33,b3:"); scanf("%f%f%f%f%f%f%f%f%f%f%f%f", &a11,&a12,&a13,&b1,&a21,&a22,&a23,&b2,&a31,&a32,&a33,&b3);if(fabs(a11)<fabs(a21)){ if(fabs(a11)>fabs(a31))A=a11;a11=a31;a31=A;B=a12;a12=a32;a32=B;C=a13;a13=a33;a33=C;D=b1;b1=b3;b3=D ;}if (fabs(a11)<fabs(a21)){if(fabs(a21)>fabs(a31)){A=a11;a11=a21;a21=A;B=a12;a12=a22;a22=B;C=a13;a13=a23;a23=C;D=b1;b1=b2;b2=D ;}elseA=a11;a11=a31;a31=A;B=a12;a12=a32;a32=B;C=a13;a13=a33;a33=C;D=b1;b1=b3;b1=D ;}printf("now a11=%f a12=%f a13=%f b1=%f\n",a11,a12,a13,b1); printf("now a21=%f a22=%f a23=%f b2=%f\n",a21,a22,a23,b2); printf("now a31=%f a32=%f a33=%f b3=%f\n",a31,a32,a33,b3);l21=a21/a11; l31=a31/a11;A22=a22-l21*a12;A23=a23-l21*a13;d1=b2-l21*b1;A32=a32-l31*a12; A33=a33-l31*a13;d2=b3-l31*b1;if(fabs(A22)>fabs(A32)){ l32=A32/A22;a=A33-l32*A23;d3=d2-l32*d1;x3=d3/a;x2=(d1-A23*x3)/A22;x1=(b1-a13*x3-a12*x2)/a11;}else l32=A22/A32;a=A23-l32*A33;d3=d1-l32*d2;x3=d3/a;x2=(d2-A33*x3)/A32;x1=(b1-a13*x3-a12*x2)/a11;printf("x1=%f x2=%f x3=%f\n",x1,x2,x3);getch(); return 0; }。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
本科实验报告
课程名称：数值计算方法B
实验项目：线性方程组的直接解法
最小二乘拟合多项式
实验地点：ZSA401
专业班级：学号：201000
学生姓名：
指导教师：李志
2012年4月13日
线性方程组的直接解法
一、实验目的和要求
实验目的：合理利用Gauss 消元法、LU 分解法或追赶法求解方程组。
实验要求：利用高斯消元法，LU 分解法或追赶法进行编程，求解题中所给的方程组。
二、实验内容和原理
实验内容：合理利用Gauss 消元法、LU 分解法或追赶法求解下列方程组：
① ⎥⎥
⎥⎦
⎤
⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡13814142210321321x x x ②⎥⎥⎥
⎥
⎦⎤
⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎦⎤⎢⎢
⎢⎢
⎢
⎣⎡--⨯-2178.4617.5911212592.1121130.6291.513
14
.59103.043
2115x x x x ③⎥⎥
⎥⎥⎥
⎥⎦⎤
⎢
⎢⎢⎢⎢⎢⎣⎡----=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-55572112112112121 n n x x x x （n=5,10,100,…）
实验原理：这个实验我选用的是高斯消元法。
高斯消元法：先按照 L ik =a ik^(k-1)/a kk^(k-1) ,
a ij^(k)=a ij^(k-1)-l ik a kj^(k-1)
[其中k=1,2,…,n-1;i=k+1,k+2,…,n;j=k+1,k+2,…,n+1]
将方程组变为上三角矩阵，再经过回代，即可求解出方程组的解。
三.计算公式
通过消元、再回代的求解方法称为高斯消元法。
特点是始终消去主对角线 下方的元素。
四、操作方法与实验步骤
#include "Stdio.h"
#define N 3 main() {
double a[N][N+1],b[N]; int i,j,k,x=0; for(i=0;i<N;i++)
{
for(j=0;j<=N;j++)
{
scanf("%lf",&a[i][j]);
}
}
for(k=0;k<N;k++)
{
for(i=x;i<N-1;i++)
{
for(j=N;j>=0;j--)
{
a[i+1][j]=a[i+1][j]-a[x][j]*a[i+1][x]/a[x][x];
}
}
x++;
}
for(i=0;i<N;i++)
{
for(j=0;j<=N;j++)
{
printf("%lf ",a[i][j]);
}
printf("\n");
}
b[2]=a[2][3]/a[2][2];
b[1]=(a[1][3]-a[1][2]*b[2])/a[1][1];
b[0]=(a[0][3]-a[0][2]*b[2]-a[0][1]*b[1])/a[0][0];
for(i=0;i<N;i++)
{
printf("x(%d)=%lf\n",i+1,b[i]);
}
getch();
}
五、实验数据记录和处理
实验结果：
六、实验结果与分析
大体来说，实现了课程设计的算法要求及功能，有很多还不能很好的处理的问题，需要我们在改进中不断完善。
做本次实验中有些地方有困难，就是获得数组中各元素的值，要用到for循环来输入各个元素的值。
还有就是将方程组化为上三角矩阵，也要用到好几个for 循环，比较容易出错。
七、讨论、心得
做实验要求我们把基础学扎实，上机实验让我又重新巩固了C语言知识。
我的实验是线性方程组的直接解法，刚开始调试代码的时候有时候就是很小的错误导致整个程序不能运行，需要我们一步一步慢慢来，经过无数次的检查程序错误的原因，以及在老师的帮助下，完成了这次实验。
最小二乘拟合多项式
一、实验内容
给定数据点（x i ,y i），用最小二乘法拟合数据的多项式，并求平方误差。
x i0 0.5 0.6 0.7 0.8 0.9 1.0 y i 1 1.75 1.96 2.19 2.44 2.71 3.00
二、计算公式
y=a0+a1x
三、结构程序设计
//最小二乘法拟合
//本程序包含两个函数，主函数与Gauss消元函数
#include<stdio.h>
#include<math.h>
#define N 10
#define M 20
double A[N][N+1];
void Gauss_eliminate(int n,double* answer){
int k,i,j,sum;
/*消元重构矩阵使之成为上三角矩阵*/
for(k=1;k<n;k++){ //k为高斯消元法上标，指示消元代数
for(i=k+1;i<=n;i++)
A[i][k]=A[i][k]/A[k][k];
for(i=k+1;i<=n;i++)
for(j=k+1;j<=n+1;j++)
A[i][j]=A[i][j]-A[k][j]*A[i][k];
}
for(i=2;i<=n;i++)
for(j=1;j<i;j++)
A[i][j]=0;
/*回代求解*/
answer[n]=A[n][n+1]/A[n][n]; //x[n]可直接解出
for(i=n-1;i>=1;i--){
sum=0;
for(j=i+1;j<=n;j++) sum+=A[i][j]*answer[j];
answer[i]=(A[i][n+1]-sum)/A[i][i];
}
}
void main(){
int k;
printf("请输入拟合次数：\n");
scanf("%d",&k);
int n=k+1;
int m,i,j;
double sum[M];
double x[M],y[M];
printf("请输入散点的个数:");
scanf("%d",&m);
printf("请创建数组x:\n");
for(i=0;i<m;i++)
scanf("%lf",&x[i]);
printf("请创建数组y:\n");
for(i=0;i<m;i++)
scanf("%lf",&y[i]);
//求系数矩阵
for(i=0;i<M;i++) sum[i]=0; //初始化
for(i=0;i<2*k+1;i++)
for(j=0;j<m;j++)
sum[i]+=pow(x[j],i);
for(i=0;i<n;i++)
for(j=0;j<m;j++)
sum[2*k+1+i]+=pow(x[j],i)*y[j];
//创建系数矩阵
int count;
for(i=1;i<=n;i++){ //先建上三角矩阵
A[i][i]=sum[2*(i-1)];
for(j=i+1,count=1;j<=n;j++){
A[i][j]=sum[2*(i-1)+count];
count++;
}
}
for(i=2;i<=n;i++) //通过对称建立下三角
for(j=1;j<i;j++)
A[i][j]=A[j][i];
for(i=1;i<=n;i++) //增广矩阵的一列
A[i][n+1]=sum[2*k+i];
//输出测试
printf("增广矩阵为:\n");
for(i=1;i<=n;i++)
{
for(j=1;j<=n+1;j++)
printf("%lf\t",A[i][j]);
printf("\n");
}
double answer[N];
Gauss_eliminate(n,answer);
/*输出解*/
for(i=1;i<=n;i++)
printf("a[%d]=%lf\t",i-1,answer[i]);
getchar();
getchar();
}
四、实验结果与讨论、心得
讨论、心得:
刚开始调试代码的时候有时候就是很小的错误导致整个程序不能运行，需要我们一步一步慢慢来，经过无数次的检查程序错误的原因，以及在老师的帮助下，完成了这次实验。
这段时间的实验课提高了我的分析问题，解决问题的能力，特别提高了对一个程序的整
体操作能力，对程序的细微之处有了明显的提高认识的地方，追求最实用的程序，弥补学习上的不足，同时认识到还应深入理解课本上的知识，学过的东西要知道理论与实践相结合，增加动手能力。