cin>>a[i][j];
}
运行结果:
下面请输入未知数的个数m=3
请按顺序输入增广矩阵a:
1234
5108
4692
4692
0-6.5-11.255.5
0-1.86265e-008-0.1153853.92308
该方程组的解为:
x[0]=-9.99999x[1]=58x[2]=-34Pressanykeytocontinue
总结:
列主元素消去法的目的是为了防止减去一个较小的数时大数淹没小数,而使结果产生较大误差,本程序关键在每次消元时找到相应列中的最大项,然后交换两行位置,在进行计算。
(2)LU分解法求解线性方程:
#include
voidsolve(floatl[][100],floatu[][100],floatb[],floatx[],intn)
{inti,j;
floatt,s1,s2;
floaty[100];
for(i=1;i<=n;i++)/*第一次回代过程开始*/
{s1=0;
for(j=1;j
{
t=-l[i][j];
s1=s1+t*y[j];
}
y[i]=(b[i]+s1)/l[i][i];}
for(i=n;i>=1;i--)/*第二次回代过程开始*/
{
s2=0;
for(j=n;j>i;j--)
{
t=-u[i][j];
s2=s2+t*x[j];
}
x[i]=(y[i]+s2)/u[i][i];
}
}
voidmain()
{floata[100][100],l[100][100],u[100][100],x[100],b[100];
inti,j,n,r,k;
floats1,s2;
for(i=1;i<=99;i++)/*将所有的数组置零,同时将L矩阵的对角值设为1*/
for(j=1;j<=99;j++)
{
l[i][j]=0,u[i][j]=0;
if(j==i)l[i][j]=1;
}
printf("inputn:
\n");/*输入方程组的个数*/
scanf("%d",&n);
printf("inputarrayA:
\n");/*读取原矩阵A*/
for(i=1;i<=n;i++)
for(j=1;j<=n;j++)
scanf("%f",&a[i][j]);
printf("inputarrayB:
\n");/*读取列矩阵B*/
for(i=1;i<=n;i++)
scanf("%f",&b[i]);
for(r=1;r<=n;r++)/*求解矩阵L和U*/
{
for(i=r;i<=n;i++)
{
s1=0;
for(k=1;k<=r-1;k++)
s1=s1+l[r][k]*u[k][i];
u[r][i]=a[r][i]-s1;
}
for(i=r+1;i<=n;i++)
{s2=0;
for(k=1;k<=r-1;k++)
s2=s2+l[i][k]*u[k][r];
l[i][r]=(a[i][r]-s2)/u[r][r];
}
}
printf("arrayL:
\n");/*输出矩阵L*/
for(i=1;i<=n;i++)
{
for(j=1;j<=n;j++)
printf("%7.3f",l[i][j]);
printf("\n");
}
printf("arrayU:
\n");/*输出矩阵U*/
for(i=1;i<=n;i++)
{
for(j=1;j<=n;j++)
printf("%7.3f",u[i][j]);
printf("\n");
}
solve(l,u,b,x,n);
printf("解为:
\n");
for(i=1;i<=n;i++)
printf("x%d=%f\n",i,x[i]);
}
运行结果:
inputn:
3
inputarrayA:
223
477
-245
inputarrayB:
31-7
arrayL:
1.0000.0000.000
2.0001.0000.000
-1.0002.0001.000
arrayU:
2.0002.0003.000
0.0003.0001.000
0.0000.0006.000
解为:
x1=2.000000
x2=-2.000000
x3=1.000000
Pressanykeytocontinue
总结:
关键是把矩阵分解为L、U两个三角矩阵,回代过程比较简单。
3:
第四章
(1)拉格朗日差值多项式;
#include
#include
#defineMAX100
voidmain()
{inti,j,k,m,n,N,mi;
floattmp,mx;
floatX[MAX][MAX],Y[MAX],x[MAX],y[MAX],a[MAX];
printf("\n输入拟合多项式的次数:
\n");
scanf("%d",&m);
printf("\n输入给定点的个数n及坐标(x,y):
\n");
scanf("%d",&N);
printf("\n");
for(i=0;iscanf("%f,%f",&x[i],&y[i]);
for(i=0;i<=m;i++)
{
for(j=i;j<=m;j++)
{
tmp=0;
for(k=0;ktmp=tmp+pow(x[k],(i+j));
X[i][j]=tmp;
X[j][i]=X[i][j];
}
}
for(i=0;i<=m;i++)
{
tmp=0;
for(k=0;ktmp=tmp+y[k]*pow(x[k],i);
Y[i]=tmp;
}
for(j=0;j{
for(i=j+1,mi=j,mx=fabs(X[j][j]);i<=m;i++)
if(fabs(X[i][j])>mx)
{
mi=i;
mx=fabs(X[i][j]);
}
if(j{
tmp=Y[j];
Y[j]=Y[mi];
Y[mi]=tmp;
for(k=j;k<=m;k++)
{
tmp=X[j][k];
X[j][k]=X[mi][k];
X[mi][k]=tmp;
}
}
for(i=j+1;i<=m;i++)
{
tmp=-X[i][j]/X[j][j];
Y[i]+=Y[j]*tmp;
for(k=j;k<=m;k++)
X[i][k]+=X[j][k]*tmp;
}
}
a[m]=Y[m]/X[m][m];
for(i=m-1;i>=0;i--)
{
a[i]=Y[i];
for(j=i+1;j<=m;j++)
a[i]-=X[i][j]*a[j];
a[i]/=X[i][i];
}
printf("\n所求的二次多项式为:
\n");
printf("P(x)=%f",a[0]);
for(i=1;i<=m;i++)
printf("+(%f)*x^%d",a[i],i);
}
运行结果:
输入拟合多项式的次数:
5
输入给定点的个数n及坐标(x,y):
3
1,2
5,3
4,2
所求的二次多项式为:
P(x)=1.980417+(0.282759)*x^1+(-0.299937)*x^2+(0.022071)*x^3+(0.016624)*x^4+(-0.0
01934)*x^5Pressanykeytocontinue
总结:
拉格朗日计算公式中,只需要知道各个点即可
4:
第五章
(1)曲线拟合:
#include
#include
#defineMAX100
voidmain()
{inti,j,k,m,n,N,mi;
floattmp,mx;
floatX[MAX][MAX],Y[MAX],x[MAX],y[MAX],a[MAX];
printf("\n输入拟合多项式的次数:
\n");
scanf("%d",&m);
printf("\n输入给定点的个数n及坐标(x,y):
\n");
scanf("%d",&N);
printf("\n");
for(i=0;iscanf("%f,%f",&x[i],&y[i]);
for(i=0;i<=m;i++)
{
for(j=i;j<=m;j++)
{
tmp=0;
for(k=0;ktmp=tmp+pow(x[k],(i+j));
X[i][j]=tmp;
X[j][i]=X[i][j];
}
}
for(i=0;i<=m;i++)
{
tmp=0;
for(k=0;ktmp=tmp+y[k]*pow(x[k],i);
Y[i]=tmp;
}
for(j=0;j{
for(i=j+1,mi=j,mx=fabs(X[j][j]);i<=m;i++)
if(fabs(X[i][j])>mx)
{
mi=i;
mx=fabs(X[i][j]);
}
if(j{
tmp=Y[j];
Y[j]=Y[mi];
Y[mi]=tmp;
for(k=j;k<=m;k++)
{
tmp=X[j][k];
X[j][k]=X[mi][k];
X[mi][k]=tmp;
}
}
for(i=j+1;i<=m;i++)
{
tmp=-X[i][j]/X[j][j];
Y[i]+=Y[j]*tmp;
for(k=j;k<=m;k++)
X[i][k]+=X[j][k]*tmp;
}
}
a[m]=Y[m]/X[m][m];
for(i=m-1;i>=0;i--)
{
a[i]=Y[i];
for(j=i+1;j<=m;j++)
a[i]-=X[i][j]*a[j];
a[i]/=X[i][i];
}
printf("\n所求的二次多项式为:
\n");
printf("P(x)=%f",a[0]);
for(i=1;i<=m;i++)
printf("+(%f)*x^%d",a[i],i);
}
输入拟合多项式的次数:
2
输入给定点的个数n及坐标(x,y):
5
1,2
5,3
2,4
8,3
-1,5
所求的二次多项式为:
P(x)=3.952280+(-0.506315)*x^1+(0.050877)*x^2Pressanykeytocontinue
5:
第六章
(1)辛普生求积方法:
#include
#defineN16/*等分数*/
floatfunc(floatx)
{floaty;
y=4.0/(1+x*x);
return(y);
}
voidgedianzhi(floaty[],floata,floath)
{inti;
for(i=0;i<=N;i++)
y[i]=func(a+i*h);
}
floatsimpson(floaty[],floath)
{floats,s1,s2;
inti;
s1=y[1];
s2=0.0;
for(i=2;i<=N-2;i=i+2)
{s1+=y[i+1];/*计算奇数项的函数值之和*/
s2+=y[i];/*计算偶数项的函数值之和*/
}
s=y[0]+y[N]+4.0*s1+2.0*s2;
return(s*h/3.0);
}
main()
{floata,b,h,s,f[N+1];
scanf("%f,%f",&a,&b);
h=(b-a)/(float)N;
gedianzhi(f,a,h);
s=simpson(f,h);
printf("s=%f\n",s);
}
运行结果:
1,3
s=1.854590
Pressanykeytocontinue
总结:
辛普生算法是一种积分方法,采用三点法插值,如果h较小的话,误差很小,因为它的插值余项
,辛普生算法比较精确,程序关键是对所取的点的取和,注意
6:
第七章
(1)改进欧拉法求解常微分方程的初值问题
#include
floatfunc(floatx,floaty)
{return(y-x);
}
floateuler(floatx0,floatxn,floaty0,intN)
{floatx,y,yp,yc,h;
inti;
x=x0;
y=y0;
h=(xn-x0)/(float)N;
for(i=1;i<=N;i++)
{yp=y+h*func(x,y);
x=x0+i*h;
yc=y+h*func(x,yp);
y=(yp+yc)/2.0;
}
return(y);
}
main()
{floatx0,xn,y0,e;
intn;
printf("\ninputn:
\n");
scanf("%d",&n);
printf("inputx0,xn:
\n");
scanf("%f,%f",&x0,&xn);
printf("inputy0:
\n");
scanf("%f",&y0);
e=euler(x0,xn,y0,n);
printf("y(%f)=%6.4f",y0,e);
}
inputn:
20
inputx0,xn:
1,6
inputy0:
2
y(2.000000)=7.0000Pressanykeytocontinue
(2)四阶龙格—库塔法
#include
floatfunc(floatx,floaty)
{return(x-y);
}
floatrunge_kutta(floatx0,floatxn,floaty0,intN)
{floatx,y,y1,y2,h,xh;
floatd1,d2,d3,d4;
inti;
x=x0;
y=y0;
h=(xn-x0)/(float)N;
for(i=1;i<=N;i++)
{xh=x+h/2;
d1=func(x,y);
d2=func(xh,y+h*d1/2.0);
d3=func(xh,y+h*d2/2.0);
d4=func(xh,y+h*d3);
y=y+h*(d1+2*d2+2*d3+d4)/6.0;
x=x0+i*h;}
return(y);
}
main()
{floatx0,xn,y0,e;
intN;
printf("\ninputn:
\n");
scanf("%d",&N);
printf("inputx0,xn:
\n");
scanf("%f,%f",&x0,&xn);
printf("inputy0:
\n");
scanf("%f",&y0);
e=runge_kutta(x0,xn,y0,N);
printf("y(%f)=%8.6f",y0,e);
}
inputn:
10
inputx0,xn:
1,2
inputy0:
5
y(5.000000)=2.833863Pressanykeytocontinue
2-2Gauss-Seidel方法
#include
#include
intgsdl(a,b,n,x,eps)
intn;
doublea[],b[],x[],eps;
{inti,j,u,v;
doublep,t,s,q;
for(i=0;i<=n-1;i++)
{u=i*n+i;
p=0.0;
x[i]=0.0;
for(j=0;j<=n-1;j++)
{v=i*n+j;
p=p+fabs(a[v]);
}
}
if(p>=fabs(a[u]))
{printf(“fail\n”);
return(-1);
}
}
p=eps+1.0;
while(p>=eps)
{for(i=0;i<=n-1;i++)
{t=x[i];s=0.0;
for(j=0;j<=n-1;j++)
{if(j!
=i)
{s=s+a[i*n+j]*x[j];
}
x[i]=(b[i]-s)/a[i*n+j];
q=fabs(x[i]-t)/(1.0+fabs(x[i]));
if(q>p)
{p=q;
}
}
}
return
(1);
}
main()
{inti;
doubleeps;
staticdoublea[4][4]={{7,2,1,-2}{9,15,3,-2}{-2,-2,11,5}{1,3,2,13}};
staticdoublex[5],b[4]={4,7