计算方法实验汇总

 1 #include<stdio.h>
 2 #include<math.h>
 3 
 4 float lmd=2;
 5 
 6 float fa(float x)
 7 {
 8     return x-lmd*(x*x*x-3*x-1)/(3*x*x-3);
 9 }
10 
11 main()
12 {
13     float x0,x1,x2,eps;
14     int count=0;
15     printf("请输入迭代初值:
");
16     scanf("%f",&x1);
17     printf("请输入精度:
");
18     scanf("%f",&eps);
19     do
20     {
21         lmd/=2;
22         x2=fa(x1);
23     }while(fabs(fa(x2))>=fabs(fa(x1)));
24     do
25     {
26         x0=x1;
27         x1=fa(x0);
28         count++;
29     }while(fabs(x1-x0)>eps);
30     printf("x*=%f
",x1);
31     printf("迭代次数:%d
",count);
32 }

牛顿下山法

 1 #include<stdio.h>
 2 #include<math.h>
 3 main()
 4 {
 5     int i,j,k=1;
 6     double sum,max;
 7     double a[3][3]={{10,-1,-2},{-1,10,-2},{-1,-1,5}};
 8     double b[3]={7.2,8.3,4.2};
 9     double x[3]={0,0,0},eps;
10     double y[3];
11     printf("请输入精度:");
12     scanf("%lf",&eps);
13     while(k<100)
14     {
15         for(i=0;i<3;i++)
16         {
17             sum=0;
18             for(j=0;j<3;j++)
19             if(i!=j)  sum+=a[i][j]*x[j];
20             y[i]=(b[i]-sum)/a[i][i];
21         }
22         max=0;
23         printf("第%d次迭代结果: ",k);
24         for(i=0;i<3;i++)
25             printf("%f ",y[i]);
26         printf("
");
27         for(j=0;j<3;j++)
28             if(max<fabs(y[j]-x[j]))  max=fabs(y[j]-x[j]);
29         if(max<eps) break;
30         else 
31         {
32             for(j=0;j<3;j++)
33                x[j]=y[j];
34             k++;
35         }
36     }
37     if(k==100)  printf("error
");
38 }

Jacobi迭代法

 1 #include<stdio.h>
 2 #include<math.h>
 3 main()
 4 {
 5     int i,j,k=1;
 6     double sum,max;
 7     double a[3][3]={{10,-1,-2},{-1,10,-2},{-1,-1,5}};
 8     double b[3]={7.2,8.3,4.2};
 9     double x[3]={0,0,0},eps;
10     double y[3]={0,0,0};
11     printf("请输入精度:");
12     scanf("%lf",&eps);
13     while(k<100)
14     {
15         for(i=0;i<3;i++)
16         {
17             sum=0;
18             for(j=0;j<3;j++)
19             if(i!=j)  sum+=a[i][j]*x[j];
20             x[i]=(b[i]-sum)/a[i][i];
21         }
22         max=0;
23         printf("第%d次迭代结果: ",k);
24         for(i=0;i<3;i++)
25             printf("%f ",x[i]);
26         printf("
");
27         for(j=0;j<3;j++)
28             if(max<fabs(y[j]-x[j]))  max=fabs(y[j]-x[j]);
29         if(max<eps) break;
30         else 
31         {
32             for(j=0;j<3;j++)
33                y[j]=x[j];
34             k++;
35         }
36     }
37     if(k==100)  printf("error
");
38 }

Gauss-Seidol迭代法

 1 #include<stdio.h>
 2 const int N=5;
 3 double x[]={0.4,0.55,0.65,0.80,0.95,1.05};
 4 double y[]={0.41075,0.57815,0.69675,0.90,1.00,1.25382};
 5 
 6 double Lagrange(double z)
 7 {
 8     double r=0;
 9     int i,j;
10     for(i=0;i<=N;i++)
11     {
12         double m=1;
13         for(j=0;j<=N;j++)
14         {
15             if(j!=i)
16                 m*=(z-x[j])/(x[i]-x[j]);
17         }
18         r+=m*y[i];
19     }
20     return r;
21 }
22 main()
23 {
24     printf("f(0.596)=%f
",Lagrange(0.596));
25     printf("f(0.99)=%f
",Lagrange(0.99));
26 }

函数插值

 1 #include<stdio.h>
 2 main()
 3 {
 4     double a[12][12],b[11],m;
 5     int i,j,k;
 6     FILE *fp;
 7     fp=fopen("D:\Gauss.txt","r");
 8     for(i=1;i<=10;i++)
 9     {
10         for(j=1;j<=10;j++)
11             fscanf(fp,"%lf",&a[i][j]);
12         fscanf(fp,"%lf",&b[i]);
13     }
14     fclose(fp);
15     for(j=1;j<=9;j++)
16     {
17         for(i=j+1;i<=10;i++)
18         {
19             m=a[i][j]/a[j][j];
20             for(k=j;k<=10;k++)
21             {
22                 a[i][k]=a[i][k]-m*a[j][k];
23             }
24             b[i]=b[i]-m*b[j];
25         }
26     }
27     for(i=10;i>0;i--)
28     {
29         double sum=0;
30         for(j=i+1;j<=10;j++)
31         {
32             sum+=a[i][j]*b[j];
33         }
34             b[i]=(b[i]-sum)/a[i][i];
35     }
36     for(i=1;i<=10;i++)
37         printf("x%d=%f
",i,b[i]);
38 }

Gauss顺序消去法

 1 #include<stdio.h>
 2 #include<math.h>
 3 main()
 4 {
 5     double a[11][11],b[11],m,temp;
 6     int i,j,k,max;
 7     FILE *fp;
 8     fp=fopen("D:\Gauss.txt","r");
 9     for(i=1;i<=10;i++)
10     {
11         for(j=1;j<=10;j++)
12             fscanf(fp,"%lf",&a[i][j]);
13         fscanf(fp,"%lf",&b[i]);
14     }
15     fclose(fp);
16     for(i=1;i<=10;i++)
17     {
18         max=i;
19         for(j=i+1;j<=10;j++)
20         {
21             if(fabs(a[j][i])>fabs(a[i][i]))
22                 max=j;
23         }
24         for(j=1;j<=10;j++)
25         {
26             temp=a[i][j];
27             a[i][j]=a[max][j];
28             a[max][j]=temp;
29         }
30         temp=b[i];
31         b[i]=b[max];
32         b[max]=temp;
33         for(k=i+1;k<=10;k++)
34         {
35             m=a[k][i]/a[i][i];
36             for(j=i;j<=10;j++)
37                 a[k][j]=a[k][j]-m*a[i][j];
38             b[k]=b[k]-m*b[i];
39         }
40     }
41     for(i=10;i>0;i--)
42     {
43         double sum=0;
44         for(j=i+1;j<=10;j++)
45             sum+=a[i][j]*b[j];
46         b[i]=(b[i]-sum)/a[i][i];
47     }
48     for(i=1;i<=10;i++)
49         printf("x%d=%f
",i,b[i]);
50 }

Gauss列主元消去法

 1 #include<stdio.h>
 2 #include<math.h>
 3 
 4 #define maxn 1000
 5 #define m 3
 6 
 7 double x[]={0,5,10,15,20,25,30,35,40,45,50,55};
 8 double y[]={0,0.000127,0.000216,0.000286,0.000344,0.000387,0.000415,0.000437,0.000451,0.000458,0.000402,0.000464};
 9 double xx[maxn],bb[maxn],gauss[maxn][maxn],temp;
10 double b[maxn];
11 
12 
13 main()
14 {
15     int i,j,k,max;
16     double mm;
17     xx[0]=12;
18     for(i=1;i<=6;i++)              
19     {
20         xx[i]=0;
21         for(k=0;k<12;k++)
22         {
23             temp=1;
24             for(j=0;j<i;j++)
25                temp*=x[k];
26             xx[i]+=temp;
27         }
28     }
29 
30     for(i=1;i<=m+1;i++)              //系数矩阵
31     {
32         b[i]=0;
33         for(k=0;k<12;k++)
34         {
35             temp=1;
36             for(j=0;j<i-1;j++)
37                 temp*=x[k];
38             b[i]+=temp*y[k];
39         }
40     }
41 
42 
43     for(i=1;i<=m+1;i++)                 //建立矩阵
44     {
45         for(j=1;j<=m+1;j++)
46             gauss[i][j]=xx[i-1+j-1];
47     }
48 
49 
50     for(i=1;i<=m+1;i++)
51     {
52         max=i;
53         for(j=i+1;j<=m+1;j++)
54         {
55             if(fabs(gauss[j][i])>fabs(gauss[i][i]))
56                 max=j;
57         }
58         for(j=1;j<=m+1;j++)
59         {
60             temp=gauss[i][j];
61             gauss[i][j]=gauss[max][j];
62             gauss[max][j]=temp;
63         }
64         temp=b[i];
65         b[i]=b[max];
66         b[max]=temp;
67         for(k=i+1;k<=m+1;k++)
68         {
69             mm=gauss[k][i]/gauss[i][i];
70             for(j=i;j<=m+1;j++)
71                 gauss[k][j]=gauss[k][j]-mm*gauss[i][j];
72             b[k]=b[k]-mm*b[i];
73         }
74     }
75     for(i=m+1;i>0;i--)
76     {
77         double sum=0;
78         for(j=i+1;j<=m+1;j++)
79             sum+=gauss[i][j]*b[j];
80         b[i]=(b[i]-sum)/gauss[i][i];
81     }
82     for(i=1;i<=m+1;i++)
83         printf("a%d=%f
",i-1,b[i]);
84 }

曲线拟合

 1 #include<stdio.h>
 2 #include<math.h>
 3 #define N 250
 4 
 5 double f(double x)
 6 {
 7     return sqrt(4-sin(x)*sin(x));
 8 }
 9 
10 void main()
11 {
12     double a=0,b=1/4;
13     double h=(b-a)/N,T=0;
14     int k;
15     for(k=1;k<N;k++)
16     {
17         T+=f(a+k*h);;
18     }
19     T=h/2*(f(a)+2*T+f(b));
20     printf("%f",T);
21 }

复合梯形公式

 1 #include<stdio.h>
 2 #include<math.h>
 3 #define N 4000
 4 
 5 double f(double x)
 6 {
 7     if(x==0)  return 1;
 8     else return sin(x)/x;
 9 }
10 
11 void main()
12 {
13     double a=0,b=1;
14     double h=(b-a)/N;
15     double s1=f(a+h/2),s2=0;
16     int k;
17     for(k=0;k<N;k++)
18     {
19         s1+=f(a+k*h+h/2);
20         s2+=f(a+k*h);
21     }
22     s1=h/6*(f(a)+4*s1+2*s2+f(b));
23     printf("%f",s1);
24 }

复合化辛普森公式

 1 #include<stdio.h>
 2 #include<math.h>
 3 
 4 double f(double x)
 5 {
 6     if(x==0)  return 1;
 7     else return sin(x)/x;
 8 }
 9 
10 void main()
11 {
12     int i,n=1,k=0;
13     double a=0,b=1,sum;
14     double T1,T2,h,S1,S2=0,C1,C2=0,R1,R2=0,eps=0.00001;
15     h=(b-a)/2;
16     T2=h*(f(a)+f(b));    
17     while(fabs(R2-R1)>eps)
18     {
19         R1=R2;
20         T1=T2;
21         S1=S2;
22         C1=C2;
23         sum=0;
24         for(i=1;i<=n;i++)
25         {
26             sum=sum+f(a+(2*i-1)*h);
27         }
28         T2=T1/2+sum*h;
29         S2=(4*T2-T1)/3;
30         C2=(16*S2-S1)/15;
31         R2=(64*C2-C1)/63;
32         n=n*2;
33         k++;
34         h=h/2;
35     }
36     printf("%f",R2);
37 }

龙贝格求积法

 1 #include<stdio.h>
 2 #include<math.h>
 3 #define N 5
 4 
 5 
 6 double f(double x,double y)
 7 {
 8     return y-2*x/y;
 9 }
10 
11 
12 main()
13 {
14     double x0=0,x1,y0=1,y1,h=0.2;
15     double k1,k2,k3,k4;
16     int n=1,i=1;
17     do
18     {
19         x1=x0+h;
20         k1=f(x0,y0);
21         k2=f(x0+h/2,y0+h/2*k1);
22         k3=f(x0+h/2,y0+h/2*k2);
23         k4=f(x1,y0+h*k3);
24         y1=y0+h/6*(k1+2*k2+2*k3+k4);
25         printf("n=%d,x%d=%f,y%d=%f
",n,i,x1,i,y1);
26         i++;
27         n++;
28         x0=x1;
29         y0=y1;
30     }while(n<N+1);
31 }

四阶龙格-库塔法