计算方法实验代码.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;
}
升级会员 