数值计算方法编程作业C语言版完整版

数值计算方法编程作业

C语言版

HEN system office room 【HEN16H-HENS2AHENS8Q8-HENH1688】

1:第二章

(1)二分法求解非线性方程：

#include<>

#include <>

#define f(x) ((x*x-1)*x-1)

void main()

{ float a,b,x,eps;

int k=0;

printf("intput eps\n");/* 容许误差 */

scanf("%f",&eps);

printf("a,b=\n");

for(;;)

{scanf("%f, %f",&a ,&b);

if(f(a)*f(b)>=0) /* 判断是否符合二分法使用的条件 */

printf("二分法不可使用,请重新输入:\n");

else break;

}

{ x=(a+b)/2;

k++;

if(f(a)*f(x)<0) /* 如果f(a)*f(x)<0，则根在区间的左半部分 */

b=x;

else if(f(a)*f(x)>0) /* 否则根在区间的右半部分 */

a=x;

else break;

}while(fabs(b-a)>eps);/*判断是否达到精度要求,若没有达到,继续循环*/

x=(a+b)/2; /* 取最后的小区间中点作为根的近似值 */

printf("\n The root is x=%f, k=%d\n",x,k);

}

运行结果：

intput eps

a,b=

2,-5

The root is x=, k=20

Press any key to continue

总结：本题关键在于两个端点的取值和误差的判断，此程序较容易。二分法收敛速度较快，但缺点是只能求解单根。

(2)牛顿法求解非线性方程：

#include <>

float f(float x) /* 定义函数f(x) */

{ return((-3*x+4)*x-5)*x+6; }

float f1(float x) /* 定义函数f(x)的导数 */

{ return (-9*x+8)*x-5; }

void main()

{ float eps,x0,x1=;

printf("input eps:\n");

scanf("%f",&eps); /* 输入容许误差 */

{ x0=x1; /* 准备下一次迭代的初值 */

x1=x0-f(x0)/f1(x0); /* 牛顿迭代 */

}while(fabs(x1-x0)>eps); /*当满足精度,输出近似根*/

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

}

程序运行结果：

总结：关键是牛顿迭代的应用，程序中最大缺点是函数及其导数已唯一给出确定不可求的随意函数的根，牛顿法比二分法收敛快，可以求重根。

2:第三章

(1)列主元素消去法求解线性方程：

#include

#define N 20

using namespace std;

void load();

float a[N][N];

int m;

int main(){

int i,j;

int c,k,n,p,r;

float x[N],l[N][N],s,d;

cout<<"下面请输入未知数的个数m=";

cin>>m;

cout<

cout<<"请按顺序输入增广矩阵a:"<

load();

for(i=0;i

{

for(j=i;j

c=(fabs(a[j][i])>fabs(a[i][i]))

j:i; /*找列最大元素*/

for(n=0;n

{s=a[i][n]; a[i][n]=a[c][n]; a[c][n]=s;} /*将列最大数防在对角线上*/

for(p=0;p

cout<

for(k=i+1;k

{

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

for(r=i;r

a[k][r]=a[k][r]-l[k][i]*a[i][r];

}

x[m-1]=a[m-1][m]/a[m-1][m-1];

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

{

d=0;

for(j=i+1;j

d=d+a[i][j]*x[j];

x[i]=(a[i][m]-d)/a[i][i]; /*求解*/

}

cout<<"该方程组的解为:"<

for(i=0;i

cout<<"x["<

//system("pause");

return 0;

}

void load()

{

int i,j;

for(i=0;i

for(j=0;j

cin>>a[i][j];

}

运行结果：

下面请输入未知数的个数m=3

请按顺序输入增广矩阵a:

1 2 3 4

5 1 0 8

4 6 9 2

该方程组的解为:

x[0]= x[1]=58 x[2]=-34 Press any key to continue

总结：列主元素消去法的目的是为了防止减去一个较小的数时大数淹没小数，而使结果产生较大误差，本程序关键在每次消元时找到相应列中的最大项，然后交换两行位置，在进行计算。

(2)LU分解法求解线性方程：

#include<>

void solve(float l[][100],float u[][100],float b[],float x[],int n)

{int i,j;

float t,s1,s2;

float y[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];

}

void main()

{float a[100][100],l[100][100],u[100][100],x[100],b[100];

int i,j,n,r,k;

float s1,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 ("input n:\n");/*输入方程组的个数*/

scanf("%d",&n);

printf ("input array A:\n");/*读取原矩阵A*/

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

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

scanf("%f",&a[i][j]);

printf ("input array B:\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("array L:\n");/*输出矩阵L*/

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

{

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

printf("% ",l[i][j]);

printf("\n");

}

printf("array U:\n");/*输出矩阵U*/

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

{

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

printf("% ",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]);

}

运行结果：

input n:

input array A:

2 2 3

4 7 7

-2 4 5

input array B:

3 1 -7

array L:

array U:

解为:

x1=

x2=

x3=

Press any key to continue

总结：关键是把矩阵分解为L、U两个三角矩阵，回代过程比较简单。3:第四章

(1)拉格朗日差值多项式;

#include<>

#define MAX 100

void main()

{ int i,j,k,m,n,N,mi;

float tmp,mx;

float X[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;i

scanf("%f,%f",&x[i],&y[i]);

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

{

for(j=i;j<=m;j++)

{

tmp=0;

for(k=0;k

tmp=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;k

tmp=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);

}

运行结果：

输入拟合多项式的次数:

输入给定点的个数n及坐标(x,y):

1,2

5,3

4,2

所求的二次多项式为:

P(x)=+*x^1+*x^2+*x^3+*x^4+

01934)*x^5Press any key to continue

总结：拉格朗日计算公式中，只需要知道各个点即可4:第五章

(1)曲线拟合：

#include<>

#define MAX 100

void main()

{ int i,j,k,m,n,N,mi;

float tmp,mx;

float X[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;i

scanf("%f,%f",&x[i],&y[i]);

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

{

for(j=i;j<=m;j++)

{

tmp=0;

for(k=0;k

tmp=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;k

tmp=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);

}

输入拟合多项式的次数:

输入给定点的个数n及坐标(x,y):

1,2

5,3

2,4

8,3

-1,5

所求的二次多项式为:

P(x)=+*x^1+*x^2Press any key to continue 5:第六章

(1)辛普生求积方法：

#include <>

#define N 16 /* 等分数 */

float func(float x)

{ float y;

y=(1+x*x);

return(y);

}

void gedianzhi(float y[],float a,float h)

{ int i;

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

y[i]=func(a+i*h);

}

float simpson(float y[],float h)

{ float s,s1,s2;

int i;

s1=y[1];

s2=;

for(i=2;i<=N-2;i=i+2)

{ s1+=y[i+1]; /* 计算奇数项的函数值之和 */

s2+=y[i]; /* 计算偶数项的函数值之和 */

}

s=y[0]+y[N]+*s1+*s2;

return(s*h/;

}

main()

{ float a,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

Press any key to continue

总结：辛普生算法是一种积分方法，采用三点法插值，如果h 较小的话，误差很小，因为它的插值余项)(2180)()4(4

ξf h a b f R ??

? ??--=，辛普生算法比较精确，程序关键是对所取的点的取和，注意 6:第七章

(1)改进欧拉法求解常微分方程的初值问题

#include <>

float func(float x,float y)

{ return(y-x);

}

float euler(float x0,float xn,float y0,int N)

{ float x,y,yp,yc,h;

int i;

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)/;

}

return(y);

}

main()

{ float x0,xn,y0,e;

int n;

printf("\ninput n:\n ");

scanf("%d",&n);

printf("input x0,xn:\n ");

scanf("%f,%f",&x0,&xn);

printf("input y0:\n ");

scanf("%f",&y0);

e=euler(x0,xn,y0,n);

printf("y(%f)=%",y0,e);

}

input n:

input x0,xn:

1,6

input y0:

y= any key to continue

(2)四阶龙格—库塔法

#include <>

float func(float x,float y)

{ return(x-y);

}

float runge_kutta(float x0,float xn,float y0,int N) { float x,y,y1,y2,h,xh;

float d1,d2,d3,d4;

int i;

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/;

d3=func(xh,y+h*d2/;

d4=func(xh,y+h*d3);

y=y+h*(d1+2*d2+2*d3+d4)/;

x=x0+i*h; }

return(y);

}

main()

{ float x0,xn,y0,e;

int N;

printf("\ninput n:\n "); scanf("%d",&N);

printf("input x0,xn:\n "); scanf("%f,%f",&x0,&xn);

printf("input y0:\n "); scanf("%f",&y0);

e=runge_kutta(x0,xn,y0,N);

printf("y(%f)=%",y0,e);

}

input n:

input x0,xn:

1,2

input y0:

y= any key to continue

2-2 Gauss-Seidel方法

#include <>

#include<>

int gsdl(a,b,n,x,eps)

int n;

double a[],b[],x[],eps;

{int i,j,u,v;

double p,t,s,q;

for(i=0;i<=n-1;i++) {u=i*n+i;

p=;

x[i]=;

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+;

while(p>=eps)

{for(i=0;i<=n-1;i++) {t=x[i];s=;

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)/+fabs(x[i]));

if(q>p)

{p=q;

}

return(1);

}

main()

{int i;

double eps;

static double a[4][4]={{7,2,1,-2}{9,15,3,-2}{-2,-2,11,5}{1,3,2,13}};

static double x[5],b[4]={4,7,-1,0};

eps=;

if(dsdl(a,b,4,x,eps)>0)

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

{printf(“x(%d)=%\n”,i,x[i]);

}