计算方法.docx
《计算方法.docx》由会员分享，可在线阅读，更多相关《计算方法.docx（16页珍藏版）》请在冰豆网上搜索。
计算方法
N-拉格朗日插值
#include
#include
#defineN100
typedefstructtag{
doublex;
doubley;
}POINT;
voidmain（）
{
inti,j,n;
doublex,temp,Ln=0;
POINTpt[N];
printf（"请输入你要输入点的个数,,1<=n<=%d：
\n",N）;
printf（"n="）;
scanf（"%d",&n）;
printf（"\n"）;
printf（"\n请输入对应的点数\n"）;
for（i=0;iscanf（"%lf,%lf",&pt[i].x,&pt[i].y）;
printf（"\n"）;
printf（"输入插值点x的值：
\n"）;
scanf（"%lf",&x）;
printf（"\n"）;
for（i=0;i{
for（j=0,temp=1;j{
if（j!
=i）
temp=temp*（x-pt[j].x）/（pt[i].x-pt[j].x）;
}
Ln=Ln+temp*pt[i].y;
}
printf（"输出：
\nLn（%lf）=%lf\n",x,Ln）;
}
分段线性插值
#include
#include
main（）
{inte,h,i,j,k,n;
floatx[10],y[10],a,sum,total,e1,e2;
printf（"inputn/n"）;
scanf（"%d",&n）;
printf（"inputX/n"）;
for（i=0;i<=n;i++）
scanf（"%f",&x[i]）;
printf（"inputY/n"）;
for（i=0;i<=n;i++）
scanf（"%f",&y[i]）;
printf（"inputx/n"）;
scanf（"%f",&a）;
for（j=1;j<=n;j++）
if（x[j-1]<=a&&a<=x[j]）
e=j;
if（e==1）
h=j;
elseif（e==n）
h=j-1;
else{e1=fabs（a-x[j-2]）;
e2=fabs（a-x[j+1]）;
if（e1h=j-1;
elseh=j;}
sum=0;
for（k=h-1;k<=h+1;k++）
{total=1;
for（j=h-1;jtotal*=（（a-x[j]）/（x[k]-x[j]））;
for（j=k+1;j<=h+1;j++）
total*=（（a-x[j]）/（x[k]-x[j]））;
sum+=y[k]*total;
}
printf（"x=%f,L=%f",a,sum）;
}
两点三次额比特插值
#include
#include
#definem4
#definen5
voidmain（）
{
inti,k;
floatx[n+1],y[n+1],yy[n+1],h,z[m];
printf（"请按行输入一系列的x值:
\n"）;
for（k=0;kscanf（"%f",&x[k]）;
printf（"请按行输入一系列的y值:
\n"）;
for（k=0;kscanf（"%f",&y[k]）;
printf（"请输入一系列的y'的值:
\n"）;
for（k=0;kscanf（"%f",&yy[k]）;
printf（"请按行输入这%d个插值点:
\n",m+1）;
for（i=0;iscanf（"%f",&z[i]）;
printf（"%f\n",z[i]）;
for（i=0;ifor（k=0;kif（z[i]>=x[k]&&z[i]<=x[k+1]）
{
h=pow（（z[i]-x[k+1]）/（x[k]-x[k+1]）,2.0）*（1+2*（z[i]-x[k]）/（x[k+1]-x[k]））*y[k]+pow（（z[i]-x[k]）/（x[k+1]-x[k]）,2.0）*（1+2*（z[i]-x[k+1]）/（x[k]-x[k+1]））*y[k+1]+pow（（z[i]-x[k+1]）/（x[k]-x[k+1]）,2.0）*（z[i]-x[k]）*yy[k]+pow（（z[i]-x[k]）/（x[k+1]-x[k]）,2.0）*（z[i]-x[k+1]）*yy[k+1];
printf（"h（%f）=%f\n",z[i],h）;
}
}
多元线性数据拟合
#include
#defineN10
main（）
{
floatx[N][N],y[N],a[N][N],b[N];
intn,i,j,m,k;
printf（"输入函数的组数n和因素的个数k:
\n"）;
scanf（"%d%d",&n,&k）;
printf（"输入一个%d*%d的二维数组（行为组数，列为因子数,最后一列为变量数）:
\n",n,k+1）;
for（i=0;i{
for（j=0;jscanf（"%f",&x[i][j]）;
scanf（"%f",&y[i]）;
}
for（i=0;ifor（j=0;j{
if（j==0&&i==0）a[i][j]=n;
elseif（i==0&&j>0）
{
a[i][j]=0;
for（m=0;ma[i][j]+=x[m][j-1];
}
elseif（i>0&&j==0）
{
a[i][j]=0;
for（m=0;ma[i][j]+=x[m][i-1];
}
elseif（i>0&&j>0）
{
a[i][j]=0;
for（m=0;ma[i][j]+=x[m][i-1]*x[m][j-1];
}
}
for（i=0;i{
if（i==0）
{
b[i]=0;
for（j=0;jb[i]+=y[j];
}
else
{
b[i]=0;
for（j=0;jb[i]+=x[j][i-1]*y[j];
}
}
printf（"数据拟合后正规方程组为:
\n"）;
for（i=0;i{
for（j=0;jprintf（"%5.3f",a[i][j]）;
printf（"%5.3f\n",b[i]）;
}
return0;
}
复化梯形公式
#include
#include
floatf（floatx）
{
returnsin（x）;
}
main（）
{
floata,b,h,f1,I;
intn,i;
printf（"请输入等分小区间个数:
\n"）;
scanf（"%d",&n）;
printf（"请输入区间:
\n"）;
scanf（"%f%f",&a,&b）;
h=（b-a）/n;
f1=0;
for（i=1;i<=n-1;i++）
{
f1+=f（a+i*h）;
}
I=h/2*（f（a）+f（b）+2*f1）;
printf（"划分%d个区间后，经复化梯形公式后结果为:
\n",2*n）;
printf（"%f\n",I）;
}
逐次分半梯形公式
#include
#include
floatf（floatx）
{
returnx/（4+x）;
}
main（）
{
floata,b,c,h,F=0,T,G;
intn=1,i;
printf（"积分值的区间:
\n"）;scanf（"%f%f",&a,&b）;
printf（"计算结果的精确度:
\n"）;scanf（"%f",&c）;
h=（b-a）/2;
F+=f（a+h）;
G=h*（f（a）+f（b））;
T=G/2+h*F;
while（fabs（T-G）>=3*c）
{
F=0;
G=T;
n=2*n;
h=h/2;
for（i=1;i<=n;i++）
F+=f（a+（2*i-1）*h）;
T=G/2+h*F;
}printf（"需要划分%d个区间，\n",n）;
printf（"x/（4+x）在区间[%f,%f]得出的积分值为:
\n",a,b）;
printf（"%f\n",T）;
return0;
}
列主元消去法
#include
#include
#defineN10
main（）
{
floata[N][N],b[N],c[N],e[N],m[N][N],f1,f2,L[N][N],U[N][N];
intd[N],n,i,j,k;
printf（"方程组的阶:
\n"）;
scanf（"%d",&n）;
printf（"方程组的系数增广矩阵,一个%d*%d的二维数组:
\n",n,n+1）;
for（i=0;i{
for（j=0;jscanf（"%f",&a[i][j]）;
scanf（"%f",&b[i]）;
}
for（k=0;k{
c[k]=0;
for（i=k;i{
if（fabs（a[i][k]）>c[k]）
{
c[k]=fabs（a[i][k]）;
d[k]=i;
}
}
if（c[k]==0）
{
printf（"系数矩阵奇异，计算停止。
\n"）;
return0;
}
if（d[k]!
=k）
{
for（j=0;j{
f1=a[d[k]][j];
a[d[k]][j]=a[k][j];
a[k][j]=f1;
}
f2=b[d[k]];
b[d[k]]=b[k];
b[k]=f2;
}
for（i=k+1;i{
m[i][k]=a[i][k]/a[k][k];
a[i][k]=m[i][k];
}
for（i=k+1;i{
for（j=k+1;ja[i][j]=a[i][j]-m[i][k]*a[k][j];
b[i]=b[i]-m[i][k]*b[k];
}
}
for（i=0;i{
for（j=0;jif（j>=i）U[i][j]=a[i][j];
elseU[i][j]=0;
}
for（i=0;i{
for（j=0;jif（jelseif（j==i）L[i][j]=1;
elseL[i][j]=0;
}
printf（"消去后方程组的增广矩阵为:
\n"）;
for（i=0;i{
for（j=0;jprintf（"%5.3f",U[i][j]）;
printf（"%5.3f\n",b[i]）;
}
for（i=n-1;i>=0;i--）
{
e[i]=0;
for（j=i+1;je[i]+=a[i][j]*b[j];
b[i]=（b[i]-e[i]）/a[i][i];
}
for（i=0;ifor（j=0;j{
if（j==i）a[i][j]=1;
elsea[i][j]=0;
}
printf（"回代求得方程组的增广矩阵为:
\n"）;
for（i=0;i{
for（j=0;jprintf（"%5.3f",a[i][j]）;
printf（"%5.3f\n",b[i]）;
}
printf（"方程组的解为:
\n"）;
for（i=0;i{
printf（"x%d=%5.3f\n",i+1,b[i]）;
}
printf（"方程组的上三角矩阵为:
\n"）;
for（i=0;i{
for（j=0;jprintf（"%5.3f",U[i][j]）;
printf（"\n"）;
}
printf（"方程组的单位下三角矩阵为:
\n"）;
for（i=0;i{
for（j=0;jprintf（"%5.3f",L[i][j]）;
printf（"\n"）;
}
return0;
}
雅克比迭代法
#include
#include
#defineN10
floatmaxn（floatx[],floaty[],intn）;
main（）
{
floata[N][N],b[N],x[N],y[N],z[N],c,sum=0;
inti,j,n,k,max;
printf（"请输入方程组的维数:
\n"）;
scanf（"%d",&n）;
printf（"请输入一个%d维方程组的增广矩阵:
\n",n）;
for（i=0;i{
for（j=0;jscanf（"%f",&a[i][j]）;
scanf（"%f",&b[i]）;
}
printf（"请输入方程组的容许误差:
\n"）;
scanf（"%f",&c）;
printf（"请输入方程组的容许最大迭代次数:
\n"）;
scanf（"%d",&max）;
printf（"请输入该%d维方程组的初始迭代向量:
\n",n）;
for（i=0;iscanf（"%f",&x[i]）;
k=0;
do
{
if（k==max）
{
printf（"输出失败，停机！
\n"）;
return0;
}
k++;
for（i=0;iy[i]=x[i];
for（i=0;i{
sum=0;
for（j=0;j{
if（j!
=i）sum+=a[i][j]*y[j];
}
x[i]=（b[i]-sum）/a[i][i];
}
}
while（maxn（x,y,n）>c）;
printf（"该%d维方程组经计算%d次后,\n得到的解为:
\n",n,k）;
for（i=0;iprintf（"x%d=%5.4f\n",i+1,x[i]）;
return0;
}
floatmaxn（floatx[],floaty[],intn）
{
floatmax=0;
inti;
for（i=0;iif（fabs（x[i]-y[i]）>max）
max=fabs（x[i]-y[i]）;
returnmax;
}
升级会员 