太原理工大学数值计算方法实验报告.doc

本科实验报告

课程名称： 计算机数值方法

实验项目： 方程求根、线性方程组的直接解
法、 线性方程组的迭代解法、代数插值和最
小二乘拟合多项式

实验地点： 行 勉 楼

专业班级： ******** 学号： *********

学生姓名： ********

指导教师： 李 誌 ， 崔 冬 华

2016年 4 月 8 日

学生姓名

实验成绩

实验名称
实验一 方程求根
实验内容和要求
熟悉使用二分法、迭代法、牛顿法、割线法等方法对给定的方程进行根的求解。选择上述方法中的两种方法求方程：f(x)=x3+4x2-10=0在[1,2]内的一个实根，且要求满足精度|x*-xn|<0.5×10-5
（1）了解非线性方程求根的常见方法，如二分法、牛顿法、割线法。
（2）加深对方程求根方法的认识，掌握算法。
（3）会进行误差分析，并能对不同方法进行比较。
实验原理
1. 二分法：如果要求已知函数 f(x) = 0 的根 (x 的解)，那先要找出一个区间 [a, b]，使得f(a)与f(b)异号。根据介值定理，这个区间内一定包含着方程式的根。求该区间的中点m=(a+b)/2，并找出 f(m) 的值。若 f(m) 与 f(a) 正负号相同，则取 [m, b] 为新的区间, 否则取 [a, m]。重复第3步和第4步，直到得到理想的精确度为止。
2. 割线法是利用牛顿迭代法的思想,在根的某个领域内，函数有直至二阶的连续导数，并且不等于0，则在领域内选取初值x0,x1,迭代均收敛。
（1） 在区间[m ,n]内输入初值x0，x1.
（2） 计算x2。x2=x1-f(x1)*(x1-x0)/(f(x1)-f(x0))
（3） x0=x1,x1=x2（4） 判断是否达到精度,若是输出x1，若否 执行（2）
主要仪器设备
HP计算机
实验记录
1.二分法
// 方程求根（二分法）.cpp : 定义控制台应用程序的入口点。
//
#include "stdafx.h"
#include"iostream"
using namespace std;
class Text
{
public:
float x, y, a, b, c, n = 0;
void Getab()
{
cout << "请输入计算区间：(以空格隔开)" << endl;
cin >> a >> b;
}
float GetY(float x)
{
y = x*x*x + 4 * x*x - 10;
return y;
}
float Calculate(float a,float b)
{
c = (a + b) / 2;
n++;
if (GetY(c) == 0 || ((b - a) / 2) < 0.000005)
{
cout << c <<"为方程的解"<< endl;
return 0;
}
if (GetY(a)*GetY(c) < 0)
{
return Calculate(a,c);
}
if (GetY(c)*GetY(b)< 0)
{
return Calculate(c,b);
}
}
};
int main()
{
cout << "方程组为：f(x)=x^3+4x^2-10=0" << endl;
float a, b;
Text text;
text.Getab();
a = text.a;
b = text.b;
text.Calculate(a, b);
return 0;
}

2.割线法：
// 方程求根（割线法）.cpp : 定义控制台应用程序的入口点。
//
#include "stdafx.h"
#include"iostream"
using namespace std;
class A
{
public:
float x0,x1,y;
float GetY(float x)
{
y= x*x*x+4*x*x-10;
return y;
}
void GetNumber()
{
cout<<"请输入两个初始近似值：(以空格隔开)" << endl;
cin >> x0;
cin >> x1;
}
void Calculate(float x0,float x1)
{
float x2;
x2 = x1 - (GetY(x1) / (GetY(x1) - GetY(x0))*(x1 - x0));
if (x2==x1)
{
cout <<x2<<"为方程的解" << endl;
}
else
{
cout << x2 << endl;
return Calculate(x1, x2);
}
}
};
int main()
{
cout << "方程组为：f(x)=x^3+4x^2-10=0" << endl;
float a, b;
A text;
text.GetNumber();
a = text.x0;
b = text.x1;
text.Calculate(a,b);
return 0;
}

心得体会
使用不同的方法，可以不同程度的求得方程的解，通过二分法计算的程序实现更加了解二分法的特点，二分法过程简单，程序容易实现，但该方法收敛比较慢一般用于求根的初始近似值，不同的方法速度不同。面对一个复杂的问题，要学会简化处理步骤，分步骤一点一点的循序处理，只有这样，才能高效的解决一个复杂问题。

实验名称
实验二 线性方程组的直接求解
实验内容和要求
合理选择利用Gauss消元法、主元素消元法、LU分解法、追赶法求解下列方程组：

（1）了解线性方程组常见的直接解法，如Guass消元法、LU分解法、追赶法。
（2）加深对线性方程组求解方法的认识，掌握算法。
（3）会进行误差分析，并能对不同方法进行比较。

实验原理
1. 高斯分解法：
⑴将原方程组化为三角形方阵的方程组：
lik=aik/akk
aij= aij- lik* akj k=1,2,…,n-1
i=k+1,k+2, …,n j=k+1,k+2, …,n+1
⑵由回代过程求得原方程组的解：
xn= ann+1/ ann
xk=( akn+1-∑akj xj)/ akk (k=n-1,n-2, …,2,1)
2. LU分解法：
将系数矩阵A转化为A=L*U, L为单位下三角矩阵，U为普通上三角矩阵，然后通过解方程组l*y=b,u*x=y,来求解x.
主要仪器设备
HP计算机
实验记录
1.高斯消元法：
#include "stdio.h"
#include "math.h"
#include <stdio.h>
double a[5][6],a0[5][6];
double l[5],tmp;
void Exchange(int i)
{
int j,l,k;
double max=a0[i][i],temp;
j=i;
for(k=i;k<=3;k++)
{
if(a0[k][i]>max)
{
max=a0[k][i];
j=k;
}
}
for(l=i;l<=4;l++)
{
temp=a0[i][l];
a0[i][l]=a0[j][l];
a0[j][l]=temp;
}
for(i=1;i<=3;i++)
{
for(j=1;j<=4;j++)
{
a[i][j]=a0[i][j];
}
}
}
void displayA()
{
int i,j;printf("\n");
for(j=1;j<=3;j++)
{
for(i=1;i<=4;i++)
printf("%lf ",a[j][i]);
printf("\n");
}
}
void main()
{
int i,j,k;
for(i=1;i<=3;i++)
{
for(j=1;j<=4;j++)
{ scanf("%lf",&a[i][j]);
a0[i][j]=a[i][j];
}
}
displayA(); printf("列主元素消元法如下");
//消元过程
k=1;
do
{
Exchange(k);
displayA();
for(i=k+1;i<=3;i++)
{
l[i]=a0[i][k]/a0[k][k];
printf("l[%i][%i]=%lf",i,k,l[i]);
for(j=k;j<=4;j++)
{
a[i][j]=a0[i][j]-l[i]*a0[k][j];
}
displayA();
}
k++;
if(k==3) break;
for(j=1;j<=3;j++)
{
for(i=1;i<=4;i++)
a0[j][i]=a[j][i];
}
}while(1);
//回代过程
l[3]=a[3][4]/a[3][3];
for(k=3;k>=1;k--)
{
tmp=0;
for(j=k+1;j<=3;j++)tmp+=a[k][j]*l[j];
l[k]=(a[k][4]-tmp)/a[k][k];
}
for(i=1;i<=3;i++)
printf("x[%i]=%lf\n",i,l[i]);
}

2.LU分解法：
#include<stdio.h>
#include<math.h>
int i,j,k,r;
double m=0,p=0;
double a[3][3];
void lu(double a[3][3])
{
for(i=1;i<=2;i++)
{
if(a[0][0]!=0)
a[i][0]=a[i][0]/a[0][0];
}
for(k=1;k<=2;k++)
{
for(j=k;j<=2;j++)
{
{
for(r=0;r<=k-1;r++)
m=m+a[k][r]*a[r][j];
} a[k][j]=a[k][j]-m;
m=0;
} for(i=k+1;i<=2;i++)
{ { for(r=0;r<=k-1;r++)
p=p+a[i][r]*a[r][k];
} a[i][k]=(a[i][k]-p)/a[k][k]; p=0; } } }
void main() {
static double a[3][3]={{1,2,3},{0,1,2},{2,4,1}};
static double b[3]={14,8,13};
double c[3];
double d[3];
double f[3][3];
double m=0;
double n=0;
int r;
int i,j;
lu(a);
printf("输出U的矩阵为\n");
for(i=0;i<=2;i++)
{
for(j=i;j<=2;j++)
{
f[i][j]=a[i][j];
printf(" %f",f[i][j]);
}
printf("\n");
}
printf("输出L的矩阵为\n");
for(i=0;i<=2;i++)
{
for(j=0;j<=i;j++)
{
if(i==j)
{
a[i][j]=1;
printf(" %f",a[i][j]);
}
else
printf(" %f",a[i][j]);
}
printf("\n");
}
c[0]=b[0];
for(i=1;i<=2;i++)
{
for(r=0;r<=i-1;r=r+1)
m=m+a[i][r]*c[r];
c[i]=b[i]-m;
m=0;
}
d[2]=c[2]/f[2][2];
for(i=1;i>=0;i=i-1)
{
for(r=2;r>i;r=r-1)
n=n+f[i][r]*d[r];
d[i]=(c[i]-n)/f[i][i];
n=0;
}
printf("所求方程组解为x1=%f, x2=%f, x3=%f",d[0],d[1],d[2]); /*根据LU分解所得两个矩阵及求解步骤计算所求X一组解*/
}

心得体会
对于求解线性方程组的各种直接方法来说各有优缺点，在所有的求解方法中都应该注意其解的精度。注意不同求解方法的不同误差求法。编写程序的时候需要一步一步慢慢来，逐步增加自己的算法知识水平和解决问题的能力。

实验名称
实验三 线性方程组的迭代求解
实验内容和要求
使用雅可比迭代法或高斯-赛德尔迭代法对下列方程组进行求解。
实验原理
雅可比迭代法：
设线性方程组
Ax=b
的系数矩阵A可逆且主对角元素a11,a22,…,ann均不为零，令
D=diag(a11,a22,…,ann)
并将A分解成
A=(A-D)+D
从而线性方程组可写成
Dx=(D-A)x+b
则有迭代公式
x(k+1)=B1x(k)+f1
其中，B1=I-D-1A,f1=D-1b。
主要仪器设备
HP计算机
实验记录
#include <stdio.h>
#include<math.h>
int main()
{
int i;
double x1[20],x2[20],x3[20];
double x10,x20,x30;
printf("please input x1,x2,x3:\n");
scanf("%lf%lf%lf",&x10,&x20,&x30);
printf(" n x1[n] x2[n] x3[n]\n");
for(i=0;i<18;i++)
{
x1[0]=x10;
x2[0]=x20;
x3[0]=x30;
x1[i+1]=0.1*x2[i]+0.2*x3[i]+0.72;
x2[i+1]=0.1*x1[i]+0.2*x3[i]+0.83;
x3[i+1]=0.2*x1[i]+0.2*x2[i]+0.84;
printf("%5d %5lf %5lf %5lf\n",i,x1[i],x2[i],x3[i]);
}
return 0;
}

心得体会
在编写算法是不熟悉，查阅了很多资料，经过反复研究和试验后实现了题目的要求，使用雅克比迭代法和高斯-赛德尔都可以得到方程的解，但相比之下，高斯-赛德尔的迭代次数要比雅克比的迭代次数少，能够更快的达到所求的解的精度。

实验名称
实验四 代数插值和最小二乘法拟合
实验内容和要求
1.学习使用拉格朗日插值法或牛顿插值法求解方法。
2.了解最小二乘法的多项式拟合的具体计算方法并且注意克服正规方程组的病态。
给定数据点（xi ,yi）如下：
xi
0
0.5
0.6
0.7
0.8
0.9
1.0
yi
1
1.75
1.96
2.19
2.44
2.71
3.00
(1) 使用拉格朗日插值法或牛顿插值法, 求f(0.856)的近似值.
(2) 用最小二乘法拟合数据的(n次)多项式，求f(0.856)的近似值.
(3) 对比、分析上两结果
实验原理
设函数在区间[a,b]上n+1互异节点x0,x1,…,xn上的函数值分别为y0,y1,…,yn,求n次插值多项式Pn(x),满足条件
Pn(xj)=yj, j=0,1,…,n
令
Ln(x)=y0l0(x)+y1l1(x)+…+ynln(x)= ∑yili(x)
其中l0(x),l1(x),…, ln(x) 为以x0,x1,…,xn为节点的n次插值基函数，则Ln(x)是一次数不超过n的多项式，且满足
Ln(xj)=yj, L=0,1,…,n
再由插值多项式的唯一性，得
Pn(x)≡Ln(x)
主要仪器设备
HP计算机
实验记录(写出实验内容中的程序代码和运行结果)(可分栏或加页)
拉格朗日插值法：
#include "stdio.h"
int main()
{

double m=1.0,a=0.856,l=0;
int i,j;
double x[6]={0.50,0.60,0.70,0.80,0.90,1.00};
double y[6]={1.75,1.96,2.19,2.44,2.71,3.00};
for(i=0;i<=5;i++)
{
for(j=0;j<=5;j++)
{
if(i==j) continue;
m=m*((a-x[j])/(x[i]-x[j]));

}
l+=y[i]*m;
m=1;
}
printf("结果为%lf",l);
return 0;
}

最小二乘法:
#include "stdio.h"
#include "math.h"
int main()
{double x[7]={0,0.5,0.6,0.7,0.8,0.9,1.0},
y[7]={1,1.75,1.96,2.19,2.44,2.71,3.00},
a0,a1,sum1=0,sum2=0,sum3=0,sum4=0,sum5=0,l,r;
int m=6,i,k;
for(i=0;i<7;i++)
{
sum1+=x[i];
sum2+=x[i]*x[i];
sum3+=y[i];
sum4+=x[i]*y[i];
sum5+=y[i]*y[i];
}
l=sum1/(m+1);
a1=(sum4-l*sum3)/(sum2-l*sum1);
a0=(sum3-sum1*a1)/(m+1);
double s=sum3*a0+sum4*a1;
r=sum5-s;
printf("y=a0+a1*x\n");
printf("a0=%f a1=%f\t\n",a0,a1,r);
double q=0.856,p;
p=a0+a1*q;
printf("y=%f\n",p);
return 0;
}

心得体会
拉格朗日插值的优点是插值多项式特别容易建立，缺点是增加节点是原有多项式不能利用，必须重新建立，即所有基函数都要重新计算，这就造成计算量的浪费。牛顿插值多项式的优点是增加节点时，原先的差商仍旧不变，仍可以使用。
数据拟合的具体作法是：对给定的数据（xi ,yi）（i=0,1,…,m），在取定的函数类中，求p(x)属于此函数类，使误差ri=p(xi)- yi (i=0,1,…,m)的平方和最小，即：
∑ri2=∑（∑p(xi)-yi）2=min
从几何意义上讲，就是寻求与给定点（xi ,yi）（i=0,1,…,m）的距离平方和为最小的曲线y=p(x)。