计算方法实验代码.docx

计算方法实验代码.docx

《计算方法实验代码.docx》由会员分享，可在线阅读，更多相关《计算方法实验代码.docx（14页珍藏版）》请在冰豆网上搜索。

计算方法实验代码

实验一拉格朗日插值法

#include

voidmain（）

{

staticfloatLx[10],Ly[10];

intn,i,j;

floatx,y=0,p;

printf（"entern="）;//输入点的个数

scanf（"%d",&n）;

printf（"enterxi\n"）;//输入点对应的x的值

for（i=1;i<=n;i++）

scanf（"%f",&Lx[i]）;

printf（"enteryi\n"）;//输入点对应的y的制

for（i=1;i<=n;i++）

scanf（"%f",&Ly[i]）;

printf（"enterx="）;//输入验证的x的值

scanf（"%f",&x）;

for（i=1;i<=n;i++）

{

p=1;

for（j=1;j<=n;j++）

{

if（i!

=j）

p=p*（x-Lx[j]）/（Lx[i]-Lx[j]）;

}

y=y+p*Ly[i];

}

printf（"y=%f\n",y）;//输入验证的x对应的y的值

}

实验二最小二乘法

#include"stdio.h"

floatgs（floata[20][20],floatb[20],intn）

{

inti,j,k,l;

floats;

k=1;

while（k!

=n+1）

{

if（a[k][k]!

=0）

{

for（i=k+1;i<=n+1;i++）

{

a[i][k]=a[i][k]/a[k][k];

b[i]=b[i]-a[i][k]*b[k];

for（j=k+1;j<=n+1;j++）

a[i][j]=a[i][j]-a[i][k]*a[k][j];

}

}

k=k+1;

}

for（k=n+1;k>=1;k--）

{

s=0;

for（l=k+1;l<=n+1;l++）

s=s+a[k][l]*b[l];

b[k]=（b[k]-s）/a[k][k];

}

return0;

}

intmain（）

{

floata[20][20]={0.0};//定义a矩阵

floatc[20][20];//定义c矩阵

floatct[20][20];//定义ct矩阵

floatx[20];//定义数组用于存放x的数据

floaty[20];//定义数组用于存放y的数据

floatb[20]={0.0};//定义b矩阵

inti,j,k,m,n;

printf（"输入所求函数的最高次数n:

\n"）;//输入n（求线性的函数输入1。

。

）

scanf（"%d",&n）;

printf（"输入测试数据的组数m:

\n"）;//输入测试数据的组数

scanf（"%d",&m）;

printf（"输入x的测试数据%d个:

\n",m）;//输入x的测试数据m个

for（i=1;i<=m;i++）

scanf（"%f",&x[i]）;

printf（"输入y的测试数据%d个:

\n",m）;//输入y的测试数据m个

for（i=1;i<=m;i++）

scanf（"%f",&y[i]）;

for（i=1;i<=m;i++）//c矩阵第一列赋值为1

c[i][1]=1.0;

//求C[][]

for（j=2;j<=n+1;j++）

for（i=1;i<=m;i++）

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

//输出C[][]

printf（"C矩阵如下：

\n"）;

for（i=1;i<=m;i++）

for（j=1;j<=n+1;j++）

{

printf（"%f",c[i][j]）;

if（j==n+1）

printf（"\n"）;

}

//求c的转置矩阵CT[][]

for（i=1;i<=m;i++）

for（j=1;j<=n+1;j++）

ct[j][i]=c[i][j];

//输出CT[][]

printf（"CT矩阵如下：

\n"）;

for（i=1;i<=n+1;i++）

for（j=1;j<=m;j++）

{

printf（"%f",ct[i][j]）;

if（j==m）

printf（"\n"）;

}

//求a[][]

for（i=1;i<=n+1;i++）

for（j=1;j<=n+1;j++）

for（k=1;k<=m;k++）

a[i][j]+=ct[i][k]*c[k][j];

//输出a[][]

printf（"a矩阵如下：

\n"）;

for（i=1;i<=n+1;i++）

for（j=1;j<=n+1;j++）

{

printf（"%f",a[i][j]）;

if（j==n+1）

printf（"\n"）;

}

//求b[]

for（i=1;i<=n+1;i++）

for（k=1;k<=m;k++）

b[i]+=ct[i][k]*y[k];

//输出b[]

printf（"b矩阵如下：

\n"）;

for（i=1;i<=n+1;i++）

printf（"%f",b[i]）;

gs（a,b,n）;//调用高斯函数求方程组的解

printf（"\n\n"）;

printf（"输出求得的函数的系数为:

\n"）;

for（i=1;i<=n+1;i++）//输出求得的函数的系数

printf（"a%d=%f",i-1,b[i]）;

printf（"\n\n"）;

return0;

}

实验三数值积分

1.定步长

#include"stdio.h"

floatf（floatx）

{

return4/（1+x*x）;

}

intmain（）

{

inti,n;

floata,b,h,xi;

floatT=0.0;

printf（"inputa,b,n:

\n"）;

scanf（"%f%f%d",&a,&b,&n）;

T=T+f（a）+f（b）;

h=（b-a）/n;

xi=a;

for（i=1;i<=n-1;i++）

{

xi=xi+h;

T=T+2*f（xi）;

}

T=T*h/2;

printf（"求得的结果为：

%f\n",T）;

return0;

}

2.变步长

#include"stdio.h"

#include"math.h"

doublef（doublex）

{

return4.0/（1+x*x）;

}

intmain（）

{

intn=1,k=1;

doublea,b,T2n,Tn;

printf（"inputa,b:

\n"）;

scanf（"%lf%lf",&a,&b）;

Tn=（f（a）+f（b））*（b-a）/2.0;

T2n=Tn/2.0+（b-a）/2.0*f（a+（b-a）/2.0）;

printf（"T1=%lf\n",Tn）;

printf（"T2=%lf\n",T2n）;

while（fabs（T2n-Tn）>=0.000001）

{

n=n*2;

Tn=T2n;

T2n=Tn/2.0;

for（k=1;k<=n;k++）

T2n=T2n+f（a+（2*k-1）*（b-a）/（2*n））*（b-a）/（2*n）;

printf（"T%d=%lf\n",2*n,T2n）;

}

printf（"\n满足条件的结果为：

%lf\n",T2n）;

return0;

}

实验四常微分方程初值问题数值解法

1.欧拉法

#include"stdio.h"

doublef（doublexi,doubleyi）{

return-yi;

}

intmain（）

{

intn;

doublea,b,y0;

printf（"inputa,b,n,y0:

\n"）;

scanf（"%lf",&a）;

scanf（"%lf",&b）;

scanf（"%d",&n）;

scanf（"%lf",&y0）;

doubleh=（b-a）/n;

doublexi=a;

doubleyi=y0;

printf（"xi\t欧拉格式yi\n"）;

printf（"%lf\t%lf\n",xi,yi）;

inti;

doublexi1,yi1;

for（i=1;i<=n;i++）

{

xi1=xi+h;

yi1=yi+h*f（xi,yi）;

xi=xi1;

yi=yi1;

printf（"%lf\t%lf\n",xi,yi）;

}

return0;

}

2.改进欧拉法

#include"stdio.h"

doublef（doublexi,doubleyi）{

return-yi;

}

intmain（）

{

intn;

doublea,b,y0;

printf（"inputa,b,n,y0:

\n"）;

scanf（"%lf",&a）;

scanf（"%lf",&b）;

scanf（"%d",&n）;

scanf（"%lf",&y0）;

doubleh=（b-a）/n;

doublexi=a;

doubleyi=y0;

printf（"xi\t改进的欧拉格式yi\n"）;

printf（"%lf\t%lf\n",xi,yi）;

inti;

doublexi1,yi1,yp,yc;

for（i=1;i<=n;i++）

{

xi1=xi+h;

yp=yi+h*f（xi,yi）;

yc=yi+h*f（xi1,yp）;

yi1=（1.0/2.0）*（yp+yc）;

xi=xi1;

yi=yi1;

printf（"%lf\t%lf\n",xi,yi）;

}

return0;

}

实验五非线性方程求解

#include"stdio.h"

#include"math.h"

doublef（doublex）

{

returnpow（x,6）-x-1;

}

doublef_（doublex）

{

return6*pow（x,5）-1;

}

intmain（）

{

doublex0,x,eps;

inti,k,N;

printf（"请输入初始x0的值：

\n"）;

scanf（"%lf",&x0）;

printf（"迭代次数：

\n"）;

scanf（"%d",&N）;

printf（"请输入误差e：

\n"）;

scanf（"%lf",&eps）;

printf（"k\txk\n"）;

printf（"0\t%.8lf\n",x0）;

x=x0-f（x0）/f_（x0）;

if（!

（f_（x0）==0））

{

for（k=1;k<=N;k++）

{

printf（"%d\t%.8lf\n",k,x）;

if（fabs（x-x0）

{

i=0;

break;

}

x0=x;

x=x0-f（x0）/f_（x0）;

}

if（k>N）

i=2;

}

else

{

i=1;

}

if（i==0）

{

printf（"\nx=%.8lf\n",x）;

printf（"f（x）=%.30lf\n",f（x））;

printf（"k=%d\n",k）;

}

elseif（i==1）

{

printf（"给定的x0导致计算中断!

!

\n"）;

}

elseif（i==2）

{

printf（"迭代N次后仍不满足精度要求!

!

\n"）;

}

return0;

}

实验六高斯消元

#include"stdio.h"

#include"math.h"

intmain（）

{

doubleA[10][10],b[10];

doublem[10];

doubleeps=0.000001;

inti,j,n,k,ik;

printf（"输入矩阵的阶数n：

\n"）;

scanf（"%d",&n）;

printf（"输入矩阵A:

\n"）;

for（i=1;i<=n;i++）

for（j=1;j<=n;j++）

scanf（"%lf",&A[i][j]）;

printf（"输入矩阵b:

\n"）;

for（i=1;i<=n;i++）

scanf（"%lf",&b[i]）;

doublemax;

doubletemp;

for（k=1;k

{

max=0;

for（i=k;i<=n;i++）

if（（max-fabs（A[i][k]））<0）

{

max=A[i][k];

ik=i;

}

if（k!

=ik）

{

for（j=k;j<=n;j++）

{

temp=A[k][j];

A[k][j]=A[ik][j];

A[ik][j]=temp;

}

temp=b[k];

b[k]=b[ik];

b[ik]=temp;

}

for（i=k+1;i<=n;i++）

{

m[i]=A[i][k]/A[k][k];

for（j=k;j<=n;j++）

{

A[i][j]=A[i][j]-A[k][j]*m[i];

}

}

for（j=k+1;j<=n;j++）

{

b[j]=b[j]-b[k]*m[j];

}

}

printf（"A上三角矩阵：

\n"）;

for（i=1;i<=n;i++）

for（j=1;j<=n;j++）

{

printf（"%12lf",A[i][j]）;

if（j==n）

printf（"\n"）;

}

printf（"b矩阵:

\n"）;

for（i=1;i<=n;i++）

{

printf（"%12lf",b[i]）;

}

b[n]=b[n]/A[n][n];

for（i=n-1;i>=1;i--）

{

for（j=n;j>i;j--）

{

b[i]=（b[i]-A[i][j]*b[j]）;

}

b[i]=b[i]/A[i][i];

}

printf（"\n求得得结果为：

\n"）;

for（i=1;i<=n;i++）

{

printf（"x%d=%lf\t",i,b[i]）;

}

return0;

}