东南大学计算方法实验报告

计算方法与实习实验报告

学院：电气工程学院

指导老师：李翠平

班级：160093

姓名：黄芃菲

学号：16009330

实习题一

实验1 拉格朗日插值法

一、方法原理

n次拉格朗日插值多项式为：L n(x)=y0l0(x)+y1l1(x)+y2l2(x)+…+y n l n(x)

n=1时，称为线性插值，L1(x)=y0(x-x1)/(x0-x1)+ y1(x-x0)/(x1-x0)=y0+(y1-x0)(x-x0)/(x1-x0)

n=2时，称为二次插值或抛物线插值，精度相对高些

L2(x)=y0(x-x1)(x-x2)/(x0-x1)/(x0-x2)+y1(x-x0)(x-x2)/(x1-x0)/(x1-x2)+y2(x-x0)(x-x1)/(x2-x0)/(x2-x1)

二、主要思路

使用线性方程组求系数构造插值公式相对复杂，可改用构造方法来插值。

对节点x i(i=0,1,…,n)中任一点x k(0<=k<=n)作一n 次多项式l k(x k)，使它在该点上取值为1，而在其余点x i(i=0,1,…,k-1,k+1,…,n)上为0，则插值多项式为L n(x)=y0l0(x)+y1l1(x)+y2l2(x)+…+y n l n(x) 上式表明：n 个点x i(i=0,1,…,k-1,k+1,…,n)都是l k(x)的零点。可求得l k

三．计算方法及过程：1.输入节点的个数n

2.输入各个节点的横纵坐标

3.输入插值点

4.调用函数，返回z

函数语句与形参说明

程序源代码如下：

形参与函数类型参数意义

int n 节点的个数

double x[n]（double *x）存放n个节点的值

double y[n]（double *y）存放n个节点相对应的函数值

double p 指定插值点的值

double fun() 函数返回一个双精度实型函数值，即插值点p

处的近似函数值

#include

#include

using namespace std;

#define N 100

double fun(double *x,double *y, int n,double p);

void main()

{int i,n;

cout<<"输入节点的个数n：";

cin>>n;

double x[N], y[N],p;

cout<<"please input xiangliang x= "<

for(i=0;i>x[i];

cout<<"please input xiangliang y= "<

for(i=0;i>y[i];

cout<<"please input LagelangrichazhiJieDian p= "<

cin>>p;

cout<<"The Answer= "<

system("pause") ;}

double fun(double x[],double y[], int n,double p)

{double z=0,s=1.0;

int k=0,i=0;

double L[N];

while(k

{ if(k==0)

{ for(i=1;i

L[0]=s*y[0];

k=k+1;}

else

{s=1.0;

for(i=0;i<=k-1;i++)s=s*((p-x[i])/(x[k]-x[i]));

for(i=k+1;i

L[k]=s*y[k];

k++;}

}

for(i=0;i

return z;

}

五．实验分析

n=2时，为一次插值，即线性插值

n=3时，为二次插值，即抛物线插值

n=1，此时只有一个节点，插值点的值就是该节点的函数值

n<1时，结果都是返回0的；这里做了n=0和n=-7两种情况

3

常用的是线性插值和抛物线插值，显然，抛物线精度相对高些

n次插值多项式Ln（x）通常是次数为n的多项式，特殊情况可能次数小于n.例如：通过三点的二次插值多项式L2(x),如果三点共线，则y=L2(x)就是一条直线，而不是抛物线，这时L2(x)是一次式。

拟合曲线光顺性差

实验2 牛顿插值法

一、方法原理及基本思路

在拉格朗日插值方法中，若增加一个节点数据，其插值的多项式需重新计算。现构造一个插值多项式Nn(x)，只需对Nn-1(x)作简单修正(如增加某项)即可得到，这样计算方便。

利用牛顿插值公式，当增加一个节点时，只需在后面多计算一项，而前面的计算仍有用；另一方面Nn(x)的各项系数恰好又是各阶差商，而各阶差商可用差商公式来计算。

由线性代数知，对任何一个不高n次的多项式P(x)=b0＋b1x＋b2x2＋…＋bnxn (幂基) ①

也可将其写成P(x)=a0＋a1(x－x0)＋a2(x－x0) (x－x1)＋…＋a n(x－x0) …(x－x n-1)

其中ai为系数，xi为给定节点，可由①求出ai 一般情况下，牛顿插值多项式Nn(x)可写成：

Nn(x)= a0＋a1(x－x0)＋a2(x－x0) (x－x1)＋…＋a n(x－x0) …(x－x n-1))

只需求出系数ai，即可得到插值多项式。

二、计算方法及过程

1.先后输入节点个数n和节点的横纵坐标，插值点的横坐标，最后输入精度e

2. 用do-while循环语句得到跳出循环时k的值

3.将k值与n-1进行比较，若在达到精度时k

函数语句与形参说明

形参与函数类型参数意义

int n 节点的个数

float x[MAX] 存放n个节点的值（从小到大）

Float y[MAX]; 存放n个节点相对应的函数值

float x0,y0 指定插值点的横纵坐标

float e 精度

程序源代码如下：

#include

#include

using namespace std;

#define MAX 100

void main()

{ float x[MAX],y[MAX];

float x0,y0,e,N1;

float N0=0;

int i,k,n;

cout<<"输入节点的个数,n=";

cin>>n;

cout<<"请输入插值节点！"<

for(i=0;i>x[i];

cout<<"请输入对应的函数值！"<

for(i=0;i>y[i];

cout<<"请输入插值点,x0=";

cin>>x0;

cout<<"请输入精度,e=";

cin>>e;

y0=1; N1=y[0]; k=0;

do { k++;

N0=N1;

y0=y0*(x0-x[k-1]);

for(i=0;i

N1=N0+y0*y[k];

} while (fabs(N1-N0) > e && k<(n-1));

if (k==(n-1)) cout<<"计算失败！";

if (k<(n-1)) cout<<"输出结果y[x0]="<

system("pause");

}

三．运行结果测试：

题一：已知f(x)=sh(x)的函数表如下：计算f(0.23)的近似值

xi 0 0.20 0.30 0.50 Sh(xi) 0 0.20134 0.30452 0.52110

实习题二

1.用牛顿法求下列方程的根：

3）lg20 x x

+-=

实验代码：

#include

#include

#define N 100

#define eps 1e-6

#define eta 1e-8

float Newton(float(*f)(float),float(*f1)(float),float x0) {

float x1,d;

int k=0;

do

{

x1=x0-(*f)(x0)/(*f1)(x0);

if(k++>N||fabs((*f1)(x1))

{

printf("\n Newton 迭代发散");

break;

}

d=fabs(x1)<1?x1-x0:(x1-x0)/x1;

x0=x1;

printf("x(%d)=%f\t",k,x0);

}

while(fabs(d)>eps&&fabs((*f)(x1))>eta);

return x1;

}

float f(float x)

{

return log10(x)+x-2;

}

float fl(float x)

{

return 1.0/(x*log(10))+1;

}

void main()

{

float x0,y0;

printf("请输入迭代初值x0\n");

scanf("%f",&x0);

printf("x(0)=%f\n",x0);

y0=Newton(f,fl,x0);

printf("方程的根为: %f\n",y0);

}

运行窗口：

实习题三

4．编写用追赶法解三对角线性方程组的程序，并解下列方程组：

2）Ax b

实验代码

#include

#include

using namespace std;

void main(){

float a=1;

float b=-4;

float c=1;

float d[11]={0,-27,-15,-15,-15,-15,-15,-15,-15,-15,-15};

float l[11];

float bb[11];

float y[11];

float x[11];

bb[1]=b;

y[1]=d[1];

int i;

for(i=2;i<11;i++){

l[i]=a/bb[i-1];

bb[i]=b-l[i]*c;

y[i]=d[i]-l[i]*y[i];

}

x[10]=y[10]/bb[10];

for(i=9;1>0;i--){

x[i]=(y[i]-c*x[i+1])/bb[i];

}

for(i=1;i<11;i++)

cout<<'x'<

}

运行窗口

#include

#include

#include

#include

#include

#define N 11

main()

{

float a[N]={0,0,1,1,1,1,1,1,1,1,1};

float b[N]={0,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4};

float c[N]={0,1,1,1,1,1,1,1,1,1,0};

float d[N]={0,-27,-15,-15,-15,-15,-15,-15,-15,-15,-15};

float x[N]={0,0,0,0,0,0,0,0,0,0,0};

float r[N]={0,0,0,0,0,0,0,0,0,0,0};

float y[N]={0,0,0,0,0,0,0,0,0,0,0};

float q;

//clrscr();

int k;

r[1]=c[1]/b[1];

y[1]=d[1]/b[1];

for(k=2;k

{

q=b[k]-r[k-1]*a[k];

r[k]=c[k]/q;

y[k]=(d[k]-y[k-1]*a[k])/q; }

y[N-1]=(d[N-1]-y[N-2]*a[N-1])/(b[N-1]-r[N-2]*a[N-1]); x[N-1]=y[N-1]; for(k=N-2;k>=1;k--) x[k]=y[k]-r[k]*x[k+1]; for(k=1;k

运行窗口

5．12310211x x x -+=- 2348311

x x x -+=- 1232106x x x -+= 123431125x x x x -+-+=

雅克比迭代法 ● 实验代码 #include #include #define eps 1e-6 #define max 100

void Jacobi(float *a,int n,float x[]) {

int i,j,k=0;

float epsilon,s;

float *y=new float[n];

for(i=0;i

while(1)

{

epsilon=0;

k++;

for(i=0;i

{

s=0;

for(j=0;j

{

if(j==i)continue;

s+=*(a+i*(n+1)+j)*x[j];

}

y[i]=(*(a+i*(n+1)+n)-s)/(*(a+i*(n+1)+i));

epsilon+=fabs(y[i]-x[i]);

}

for(i=0;i

if(epsilon

{printf("迭代次数为:%d\n",k);return;}

if(k>=max)

{printf("迭代发散");return;}

}

delete y;

}

void main()

{

int i;

float a[4][5]={10,-1,2,0,-11,0,8,-1,3,-11,2,-1,10,0,6,-1,3,-1,11,25};

float x[4];

Jacobi(a[0],4,x);

for(i=0;i<4;i++)printf("x[%d]=%f\n",i,x[i]);

}

运行窗口

高斯-塞德尔迭代法

实验代码

#include

#include

#define N 500

void main()

{

int i;

float x[4];

float c[4][5]={10,-1,2,0,-11,0,8,-1,3,-11,2,-1,10,0,6,-1,3,-1,11,25};

void GaussSeidel(float *,int,float[]);

GaussSeidel(c[0],4,x);

for(i=0;i<4;i++)printf("x[%d]=%f\n",i,x[i]);

}

void GaussSeidel(float *a,int n,float x[])

{

int i,j,k=1;

float d,dx,eps;

for(i=0;i

while(1)

{

eps=0;

for(i=0;i

{

d=0;

for(j=0;j

{

if(j==i)continue;

d+=*(a+i*(n+1)+j)*x[j];

}

dx=(*(a+i*(n+1)+n)-d)/(*(a+i*(n+1)+i));

eps+=fabs(dx-x[i]);

x[i]=dx;

}

if(eps<1e-6){printf("迭代次数为:%d\n",k);return;}

if(k>N)

{

printf("迭代发散\n");

return;

}

k++;

}

}

●运行窗口

实习题四

2. 按下列数据

Xi 0.30 0.42 0.50 0.58 0.66 0.72

Yi 1.04403 1.08462 1.11803 1.15603 1.19817 1.23223 作5次插值，并求X1=0.46，X2=0.55，X3=0.60时的函数近似值。

●实验代码

#include

#define N 5

void Difference(float x[],float y[],int n)

{

float *f=new float[n+1];

int k,i;

for(k=1;k<=n;k++)

{

f[0]=y[k];

for(i=0;i

f[i+1]=(f[i]-y[i])/(x[k]-x[i]);

y[k]=f[k];

}

delete f;

return;

}

void main()

{

int i;

float a,b,c,varx=0.46,vary=0.55,varz=0.60;

float x[N+1]={0.30,0.42,0.50,0.58,0.66,0.72};

float y[N+1]={1.04403,1.08462,1.11803,1.15603,1.19817,1.23223};

Difference(x,y,N);

a=y[N];b=y[N];c=y[N];

for(i=N-1;i>=0;i--)

a=a*(varx-x[i])+y[i];

for(i=N-1;i>=0;i--)

b=b*(vary-x[i])+y[i];

for(i=N-1;i>=0;i--)

c=c*(varz-x[i])+y[i];

printf("Nn(%f)=%f\n",varx,a);

printf("Nn(%f)=%f\n",vary,b);

printf("Nn(%f)=%f\n",varz,c);

}

运行窗口

实习题5

1.

1）试分别用抛物线y=a+bx+cxx和指数曲线y=ae(bx)拟合下列数据：

Xi 1 1.5 2 2.5 3 3.5

Yi 33.4 79.50 122.65 159.05 189.15 214.15 Xi 4 4.5 5 5.5 6 6.5

Yi 238.65 252.50 267.55 280.50 296.65 301.40

XI 7 7.5 8

YI 310.40 318.15 325.15 实验运行代码：

//最小二乘法拟合曲线

#include

#include

using namespace std;

const int n=15;//数据对个数

const int m=3;//线性无关函数个数

class ColPivot

{

friend class Approx;

private:

double c[m][m+1];

double x[m];

public:

ColPivot()

{

for(int i=0;i

{

for(int j=0;j

x[i]=0;

}

}

void Cal();

};

void ColPivot::Cal()

{

int i,j,t,k;

double p;

for(i=0;i

{

k=i;

for(j=i+1;j

{

if(fabs(c[j][i])>fabs(c[k][i])) k=j;

}

if(k!=j)

{

for(j=i;j

{

p=c[i][j];

c[i][j]=c[k][j];

c[k][j]=p;

}

}

for(j=i+1;j

{

p=c[j][i]/c[i][i];

for(t=i;t

{

c[j][t]-=p*c[i][t];

}

}

}

for(i=m-1;i>=0;i--)

{

p=c[i][m];

for(j=m-1;j>i;j--)

{

p-=c[i][j]*x[j];

}

if(c[i][i]==0) x[i]=0;

else x[i]=p/c[i][i];

}

}

class Approx

{

private:

double x[n],y[n],a[3],b[2];

public:

void GetValue(double[],double[]);

void Para();

void Expo();

void Show();

};

void Approx::GetValue(double a[n],double b[n]) {

for(int i=0;i

{

x[i]=a[i];

y[i]=b[i];

}

}

void Approx::Para()

{

int i,j,t;

ColPivot col;

for(i=0;i<3;i++)

{

for(j=0;j<3;j++)

{

for(t=0;t

}

for(t=0;t

}

col.Cal();

for(i=0;i<3;i++) a[i]=col.x[i];

}

void Approx::Expo()

{

int i,j,t;

ColPivot col;

for(i=0;i<2;i++)

{

for(j=0;j<2;j++)

{

for(t=0;t

}

for(t=0;t

}

col.Cal();

b[0]=exp(col.x[0]);

b[1]=col.x[1];

}

void Approx::Show()

{

cout<<"用抛物线拟合："<

cout<<"a="<

cout<<"b="<

cout<<"c="<

cout<<"用指数曲线拟合："<

cout<<"a="<

cout<<"b="<

}

int main()

{

Approx ap;

double x[15];

for(int i=0;i<15;i++) x[i]=1+i*0.5;

double

y[15]={33.4,79.5,122.65,159.05,189.15,214.15,238.65,252.5,267.55,280.5,296.65,301.4,310. 4,318.15,325.15};

ap.GetValue(x,y);

ap.Para();

ap.Expo();

ap.Show();

return 0;

}