数值分析常见算法C++实现

1.1-有效数字丢失现象观察

 1 #include<bits./stdc++.h>
 2 using namespace std;
 3 double f1(double x)
 4 {
 5     return sqrt(x)*(sqrt(x+1)-sqrt(x));
 6 }
 7 double f2(double x)
 8 {
 9     return sqrt(x)/(sqrt(x+1)+sqrt(x));
10 }
11 int main() 
12 {
13     double t=1.0;
14     while (t<=1e14)
15     {
16         printf("At t=%.lf  %.19lf %.19lf\n",t,f1(t),f2(t));
17         t*=10.0;
18     }
19     printf("sqrt(x+1) = %.19lf  sqrt(x) = %.19lf\n",sqrt(t+1),sqrt(t));
20     printf("diff = %.19lf",sqrt(t+1)-sqrt(t));
21     printf("sum =  %.19lf",sqrt(t+1)+sqrt(t)); 
22 }
View Code

1.2-n次插值的Lagrange形式

 1 #include<bits./stdc++.h>
 2 using namespace std;
 3 int n;
 4 double x[1000],y[1000],xx,yy;
 5 int main()
 6 {
 7     scanf("%d",&n);
 8     for (int i=0;i<=n;i++) scanf("%lf %lf",&x[i],&y[i]);
 9     scanf("%lf",&xx);
10     for (int i=0;i<=n;i++)
11     {
12         double t=1.0;
13         for (int k=0;k<=n;k++)
14             if (i!=k) t=t*(xx-x[k])/(x[i]-x[k]);
15         yy=yy+t*y[i];
16     }
17     printf("%.19lf",yy);
18 }
View Code

1.3-n次插值Newton形式

 1 #include<bits./stdc++.h>
 2 using namespace std;
 3 int n;
 4 double x[1000],y[1000],v[1000],xx,yy,w;
 5 int main()
 6 {
 7     int n;
 8     scanf("%d",&n);
 9     for (int i=0;i<=n;i++) scanf("%lf%lf",&x[i],&y[i]);
10     scanf("%lf",&xx);
11     w=1.0,v[0]=y[0],yy=y[0];
12     for (int k=1;k<=n;k++)
13     {
14         v[k]=y[k];
15         for (int i=0;i<=k-1;i++)
16             v[k]=(v[i]-v[k])/(x[i]-x[k]);
17         w=w*(xx-x[k-1]);
18         yy=yy+w*v[k];
19     }
20     printf("%.19lf",yy);
21 }
View Code

1.4-三次样条插值

 1 #include<bits./stdc++.h>
 2 using namespace std;
 3 int n,type;
 4 double x[1000],y[1000],h[1000],xx,yy,w;
 5 double sqr(double x)
 6 {
 7     return x*x;
 8 } 
 9 double m[1000],A[1000],B[1000],a[1000],b[1000];
10 int main()
11 {
12     int n;
13     scanf("%d",&n);
14     for (int i=0;i<=n;i++) scanf("%lf%lf",&x[i],&y[i]);
15     scanf("%lf",&xx); 
16     scanf("%d",&type); 
17     for (int i=0;i<=n-1;i++) h[i]=x[i+1]-x[i];
18     if (type==1)
19     {
20         scanf("%lf%lf",&m[0],&m[n]);
21         a[0]=0.0; b[0]=2.0*m[0];
22         a[n]=1.0; b[n]=2.0*m[n];
23     }else if (type==2)
24     {
25         a[0]=1.0;b[0]=3.0/ h[0]   *(y[1] - y[0]);
26         a[n]=0.0;b[n]=3.0/(h[n-1])*(y[n] - y[n-1]);
27     }
28     for (int i=1;i<=n-1;i++)
29     {
30         a[i]=h[i-1]/(h[i-1]+h[i]);
31         b[i]=3.0*((1-a[i]) / h[i-1] * (y[i]   - y[i-1]) + 
32               a[i]     / h[i]   * (y[i+1] - y[i])); 
33     }
34     A[0]=-0.5*a[0];
35     B[0]=0.5*b[0];
36     for (int i=1;i<=n;i++)
37     {
38         A[i]=-a[i]/(2 + (1 - a[i]) * A[i-1]);
39         B[i]=(b[i]-(1-a[i])*B[i-1])/(2+(1-a[i])* A[i-1]);
40     }
41     m[n+1]=0;
42     for (int i=n;i>=0;i--) m[i]=A[i]*m[i+1]+B[i]; 
43     for (int i=0;i<=n-1;i++)
44         if (xx>=x[i] && xx<=x[i+1])
45         {
46             yy=( (1 + 2 * (xx-x[i]) / (x[i+1]-x[i]))  * sqr((xx-x[i+1]) / (x[i]-x[i+1])) * y[i]    )+
47              ( (1 + 2 * (xx-x[i+1])/ (x[i]-x[i+1])) * sqr((xx-x[i])   / (x[i+1]-x[i])) * y[i+1]  )+
48             ( (xx-x[i]  ) * sqr( (xx-x[i+1]) / (x[i]-x[i+1]) ) * m[i]   )+
49              ( (xx-x[i+1]) * sqr( (xx-x[i])   / (x[i+1]-x[i]) ) *m[i+1]  );
50             printf("%.19lf\n",yy);
51         }
52 }
View Code

1.5-单变量数据拟合（最小二乘法）

 1 #include<bits./stdc++.h>
 2 using namespace std;
 3 int n,type;
 4 double x[1000],y[1000],sx,sy,sxx,sxy,a,b;
 5 int main()
 6 {
 7     scanf("%d",&n);
 8     for (int i=1;i<=n;i++) scanf("%lf%lf",&x[i],&y[i]);
 9     for (int i=1;i<=n;i++)
10     {
11         sx+=x[i];
12         sy+=y[i];
13         sxx+=x[i]*x[i];
14         sxy+=x[i]*y[i];
15     }
16     double nn=n;
17     a=(sxx*sy-sx*sxy)/(nn*sxx-sx*sx);
18     b=(nn*sxy-sx*sy)/(nn*sxx-sx*sx);
19     printf("%.19lf %.19lf\n",a,b);
20 }
View Code

1.6-自动选取步长梯形法

 1 #include<bits./stdc++.h>
 2 using namespace std;
 3 double a,b,c,h,t1,t0,s,n;
 4 double f(double x)
 5 {
 6     return 2/(1+x*x);
 7 }
 8 int main()
 9 {
10     scanf("%lf%lf%lf",&a,&b,&c);
11     h=(b-a)/2.0;
12     t1=(f(a)+f(b))*h;
13     n=1;
14     while (1)
15     {
16         t0=t1;
17         s=0;
18         for (int k=1;k<=n;k++) s=s+f(a+(2*k-1)*h/n);
19         t1=t0/2.0+s*h/n;
20         if (abs(t1-t0)<3*c)
21         {
22             printf("%.19lf",t1);
23             break;
24         }else n=2.0*n;
25     }
26 }
View Code

1.7-Romberg 求积法

 1 #include<bits./stdc++.h>
 2 using namespace std;
 3 double a,b,c;
 4 double t[1000][1000];
 5 int k;
 6 double f(double x)
 7 {
 8     return sqrt(x);
 9 }
10 double why(double p,double q,int k)
11 {
12     double sum=0;
13     double kk=k;
14     for (double i=p;i<=q;i+=1.0)
15     {
16         sum+=f(a+(2*i-1)*(b-a)/pow(2,kk));
17     } 
18     return sum;
19 }
20 int main()
21 {
22     scanf("%lf%lf%lf",&a,&b,&c);
23     t[0][0]=(b-a)/2.0*(f(a)+f(b));
24     k=1;
25     while (1)
26     {
27     
28         t[0][k]=0.5*(t[0][k-1]+(b-a)/(pow(2,k-1))*why(1,pow(2,k-1),k) );
29         for (int m=1;m<=k;m++)
30         {
31             t[m][k-m]=(pow(4,m)*t[m-1][k-m+1]-t[m-1][k-m])/(pow(4,m)-1);
32         }
33         if (abs(t[k][0]-t[k-1][0])<c) 
34         {
35             printf("%.19lf",t[k][0]);
36             break;
37         }else k+=1.0;
38     }
39 }
View Code

2.1-顺序高斯消去法

 1 #include<bits./stdc++.h>
 2 using namespace std;
 3 int n,m;
 4 double a[100][100],b[100],c,x[100],t;
 5 double sum(int i)
 6 {
 7     double su=0;
 8     for (int j=i+1;j<=n;j++) su+=a[i][j]*x[j];
 9     return su;
10 }
11 int main()
12 {
13     scanf("%d%d",&n,&m);
14     for (int i=1;i<=n;i++)
15     for (int j=1;j<=m;j++) scanf("%lf",&a[i][j]);
16     for (int i=1;i<=n;i++) scanf("%lf",&b[i]);
17     scanf("%lf",&c); 
18     for (int k=1;k<=n-1;k++)
19     {
20         if (abs(a[k][k])<=c) 
21         {
22             printf("求解失败\n");
23             return 0;
24         }
25         for (int i=k+1;i<=n;i++)
26         {
27             t=a[i][k]/a[k][k];
28             b[i]=b[i]-t*b[k];
29             for (int j=k+1;j<=n;j++) a[i][j]=a[i][j]-t*a[k][j];
30         }
31     }
32     if (abs(a[n][n])<=c)
33     {
34         printf("求解失败\n");
35         return 0;
36     }else x[n]=b[n]/a[n][n];
37     for (int i=n-1;i>=1;i--) x[i]=(b[i]-sum(i))/a[i][i];
38     for (int i=1;i<=n;i++) printf("%.19lf ",x[i]);
39     cout<<endl;
40 }
View Code

2.2-列主元高斯消去法

 1 #include<bits./stdc++.h>
 2 using namespace std;
 3 int n,m,ik;
 4 double a[100][100],b[100],c,x[100],t;
 5 double sum(int i)
 6 {
 7     double su=0;
 8     for (int j=i+1;j<=n;j++) su+=a[i][j]*x[j];
 9     return su;
10 }
11 int main()
12 {
13     scanf("%d%d",&n,&m);
14     for (int i=1;i<=n;i++)
15     for (int j=1;j<=m;j++) scanf("%lf",&a[i][j]);
16     for (int i=1;i<=n;i++) scanf("%lf",&b[i]);
17     scanf("%lf",&c); 
18     for (int k=1;k<=n-1;k++)
19     {
20         t=-1e9;
21         for (int i=k;i<=n;i++) 
22             if (abs(a[i][k])>t)
23             {
24                 t=abs(a[i][k]);
25                 ik=i;
26             }
27         if (t<=c) {
28             printf("求解失败\n");
29             return 0;
30         }
31         if (ik!=k)
32         {
33             for (int i=1;i<=m;i++) swap(a[ik][i],a[k][i]);
34             swap(b[ik],b[k]);
35         }
36         for (int i=k+1;i<=n;i++)
37         {
38             t=a[i][k]/a[k][k];
39             b[i]=b[i]-t*b[k];
40             for (int j=k+1;j<=n;j++) a[i][j]=a[i][j]-t*a[k][j];
41         }
42     }
43     if (abs(a[n][n])<=c)
44     {
45         printf("求解失败\n");
46         return 0;
47     }else x[n]=b[n]/a[n][n];
48     for (int i=n-1;i>=1;i--) x[i]=(b[i]-sum(i))/a[i][i];
49     for (int i=1;i<=n;i++) printf("%.19lf ",x[i]);
50     cout<<endl;
51 }
View Code

2.3-全主元高斯消去法

 1 #include<bits./stdc++.h>
 2 using namespace std;
 3 int n,m,ik,jk,d[100];
 4 double a[100][100],b[100],c,x[100],t,z[100];
 5 double sum(int i)
 6 {
 7     double su=0;
 8     for (int j=i+1;j<=n;j++) su+=a[i][j]*z[j];
 9     return su;
10 }
11 int main()
12 {
13     scanf("%d%d",&n,&m);
14     for (int i=1;i<=n;i++)
15     for (int j=1;j<=m;j++) scanf("%lf",&a[i][j]);
16     for (int i=1;i<=n;i++) scanf("%lf",&b[i]);
17     scanf("%lf",&c); 
18     for (int i=1;i<=n;i++) d[i]=i;
19     for (int k=1;k<=n-1;k++)
20     {
21         t=-1e9;
22         for (int i=k;i<=n;i++) 
23         for (int j=k;j<=n;j++)
24         {
25             if (abs(a[i][j])>t)
26             {
27                 t=abs(a[i][j]);
28                 ik=i;
29                 jk=j;
30             }
31         }
32         if (t<=c) 
33         {
34             printf("求解失败\n");
35             return 0;
36         }
37         if (ik!=k)
38         {
39             for (int i=1;i<=m;i++) swap(a[ik][i],a[k][i]);
40             swap(b[ik],b[k]);
41         }
42         if (jk!=k)
43         {
44             for (int i=1;i<=n;i++) swap(a[i][jk],a[i][k]);
45             swap(d[jk],d[k]);
46         }
47         for (int i=k+1;i<=n;i++)
48         {
49             t=a[i][k]/a[k][k];
50             b[i]=b[i]-t*b[k];
51             for (int j=k+1;j<=n;j++) a[i][j]=a[i][j]-t*a[k][j];
52         }
53     }
54     z[n]=b[n]/a[n][n];
55     for (int i=n-1;i>=1;i--) z[i]=(b[i]-sum(i))/a[i][i];
56     for (int j=1;j<=n;j++) x[d[j]]=z[j];
57     for (int i=1;i<=n;i++) printf("%.19lf ",x[i]);
58     cout<<endl;
59 }
View Code

2.4-LU 直接分解法

 1 #include<bits./stdc++.h>
 2 using namespace std;
 3 int n,m,ik,jk,d[100];
 4 double a[100][100],b[100],c,x[100],y[100];
 5 double l[100][100],u[100][100];
 6 double sum2(int p,int i)
 7 {
 8     double su=0;
 9     for (int j=1;j<=p;j++) su+=l[i][j]*y[j];
10     return su;
11 }
12 double sum3(int p,int i)
13 {
14     double su=0;
15     for (int j=p;j<=n;j++) su+=u[i][j]*x[j];
16     return su;
17 }
18 int main()
19 {
20     scanf("%d%d",&n,&m);
21     for (int i=1;i<=n;i++)
22     for (int j=1;j<=m;j++) scanf("%lf",&a[i][j]);
23     for (int i=1;i<=n;i++) scanf("%lf",&b[i]);
24     scanf("%lf",&c);     
25     for(int k=1;k<=n;k++)
26     {
27         for(int j=k;j<=n;j++)
28         {      
29             double su=0;
30             for(int r=1;r<=k-1;r++)  su+=(l[k][r]*u[r][j]);
31             u[k][j]=a[k][j]-su;
32         }
33         for(int i=k;i<= n;i++)
34         {  double su=0;
35             for(int r=1;r<=k-1;r++)  su+=(l[i][r]*u[r][k]);
36             l[i][k]=(a[i][k]-su)/u[k][k];
37         }
38     }
39     y[1]=b[1];
40     for (int i=2;i<=n;i++) y[i]=b[i]-sum2(i-1,i);
41     x[n]=y[n]/u[n][n];
42     for (int i=n-1;i>=1;i--) x[i]=(y[i]-sum3(i+1,i))/u[i][i];
43     for (int i=1;i<=n;i++) printf("%.19lf\n",x[i]);
44     cout<<endl;
45 }
View Code

2.5-Jacobi 迭代算法

 1 #include<bits./stdc++.h>
 2 using namespace std;
 3 int n,k,M;
 4 double a[100][100],b[100],c,x[100],y[100],g[100];
 5 int main()
 6 {
 7     scanf("%d",&n);
 8     for (int i=1;i<=n;i++)
 9     { 
10         for (int j=1;j<=n;j++) scanf("%lf",&a[i][j]);
11         scanf("%lf",&b[i]);
12     } 
13     for (int i=1;i<=n;i++) scanf("%lf",&y[i]);
14     scanf("%lf%d",&c,&M);
15     k=1;
16     for (int i=1;i<=n;i++)
17     {
18         if (abs(a[i][i])<c)
19         {
20             printf("求解失败2\n");
21             return 0;
22         }
23         double t=a[i][i];
24         for (int j=1;j<=n;j++)
25         {
26             a[i][j]=-a[i][j]/t;
27         } 
28         a[i][i]=0;
29         g[i]=b[i]/t;
30     }
31     while (1)
32     {     
33         for (int i=1;i<=n;i++)
34         {
35             double sum=0;
36             for (int j=1;j<=n;j++)  if (i!=j) sum+=a[i][j]*y[j];
37             x[i]=g[i]+sum;
38         }
39         double sum=0;
40         for (int i=1;i<=n;i++) 
41         {
42             sum+=abs(x[i]-y[i]);
43         }
44         if (sum<c)
45         {
46             for (int i=1;i<=n;i++) printf("%.19lf\n",x[i]);
47             printf("%d\n",k);
48             return 0;
49         }else 
50         if (k<M)
51         {
52             k=k+1;
53             for (int i=1;i<=n;i++) y[i]=x[i];
54             for (int i=1;i<=n;i++) printf("%.15lf  ",y[i]);
55             printf("\n");
56         }else
57         {
58             printf("求解失败1\n");
59             return 0;
60         }
61     }
62 }
View Code

2.6-Seidel 迭代算法

 1 #include<bits./stdc++.h>
 2 using namespace std;
 3 int n,k,M;
 4 double a[100][100],b[100],c,x[100],y[100],g[100];
 5 int main()
 6 {
 7     scanf("%d",&n);
 8     for (int i=1;i<=n;i++)
 9     { 
10         for (int j=1;j<=n;j++) scanf("%lf",&a[i][j]);
11         scanf("%lf",&b[i]);
12     } 
13     for (int i=1;i<=n;i++) scanf("%lf",&y[i]);
14     scanf("%lf%d",&c,&M);
15     k=1;
16     for (int i=1;i<=n;i++) x[i]=y[i]; 
17     
18     for (int i=1;i<=n;i++)
19     {
20         if (abs(a[i][i])<c)
21         {
22             printf("求解失败2\n");
23             return 0;
24         }
25         double t=a[i][i];
26         for (int j=1;j<=n;j++)
27         {
28             a[i][j]=-a[i][j]/t;
29         } 
30         a[i][i]=0;
31         g[i]=b[i]/t;
32     }
33     while (1)
34     {     
35         for (int i=1;i<=n;i++)
36         {
37             double sum=0;
38             for (int j=1;j<=n;j++)  if (i!=j) sum+=a[i][j]*x[j];
39             x[i]=g[i]+sum;
40         }
41         double sum=0;
42         for (int i=1;i<=n;i++) 
43         {
44             sum+=abs(x[i]-y[i]);
45         }
46         if (sum<c)
47         {
48             for (int i=1;i<=n;i++) printf("%.19lf\n",x[i]);
49             printf("%d\n",k);
50             return 0;
51         }else 
52         if (k<M)
53         {
54             k=k+1;
55             for (int i=1;i<=n;i++) y[i]=x[i];
56             for (int i=1;i<=n;i++) printf("%.15lf  ",y[i]);
57             printf("\n");
58         }else
59         {
60             printf("求解失败1\n");
61             return 0;
62         }
63     }
64 }
View Code

2.7-松弛迭代算法

 1 #include<bits./stdc++.h>
 2 using namespace std;
 3 int n,k,M;
 4 double a[100][100],b[100],c,x[100],y[100],g[100],w;
 5 int main()
 6 {
 7     scanf("%d",&n);
 8     for (int i=1;i<=n;i++)
 9     { 
10         for (int j=1;j<=n;j++) scanf("%lf",&a[i][j]);
11         scanf("%lf",&b[i]);
12     } 
13     for (int i=1;i<=n;i++) scanf("%lf",&y[i]);
14     scanf("%lf",&w);
15     scanf("%lf%d",&c,&M);
16     k=1;
17     for (int i=1;i<=n;i++) x[i]=y[i]; 
18     
19     for (int i=1;i<=n;i++)
20     {
21         if (abs(a[i][i])<c)
22         {
23             printf("求解失败2\n");
24             return 0;
25         }
26         double t=a[i][i];
27         for (int j=1;j<=n;j++)
28         {
29             a[i][j]=-w*a[i][j]/t;
30         } 
31         a[i][i]=1-w;
32         g[i]=w*b[i]/t;
33     }
34     while (1)
35     {     
36         for (int i=1;i<=n;i++)
37         {
38             double sum=0;
39             for (int j=1;j<=n;j++)  if (i!=j) sum+=a[i][j]*x[j];
40             x[i]=g[i]+sum;
41         }
42         double sum=0;
43         for (int i=1;i<=n;i++) 
44         {
45             sum+=abs(x[i]-y[i]);
46         }
47         if (sum<c)
48         {
49             for (int i=1;i<=n;i++) printf("%.19lf\n",x[i]);
50             printf("%d\n",k);
51             return 0;
52         }else 
53         if (k<M)
54         {
55             k=k+1;
56             for (int i=1;i<=n;i++) y[i]=x[i];
57             for (int i=1;i<=n;i++) printf("%.15lf  ",y[i]);
58             printf("\n");
59         }else
60         {
61             printf("求解失败1\n");
62             return 0;
63         }
64     }
65 }
View Code

2.8-对分法求根

 1 #include<bits./stdc++.h>
 2 using namespace std;
 3 double d,c,a[1000],b[1000],A,B,x[1000];
 4 double f(double x)
 5 {
 6     return x*x*x-3*x+1.0; 
 7 } 
 8 int k;
 9 double ff;
10 int main()
11 {
12     scanf("%lf%lf",&d,&c);
13     scanf("%lf%lf",&A,&B);
14     a[0]=A;b[0]=B;k=0;
15     while (1)
16     {
17         x[k]=(a[k]+b[k])/2;
18         ff=f(x[k]);
19         if (abs(ff)<d) 
20         {
21             printf("%.19lf\n",x[k]);
22             break;
23         }
24         if (f(a[k])*f(x[k]) < 0) 
25         {
26             a[k+1]=a[k];
27             b[k+1]=x[k];
28         }else
29         {
30             a[k+1]=x[k];
31             b[k+1]=b[k];
32         }
33         if (b[k+1]-a[k+1]<=c)
34         {
35             printf("%.19lf\n",(a[k+1]+b[k+1])/2);
36         }else k=k+1;
37     }
38 }
View Code

2.9-松弛法迭代求根

#include<bits./stdc++.h>
using namespace std;
double d,c,A,B,x[1000],w[1000];
double f(double x)
{
    return (x*x*x+1.0)/3.0; 
} 
double ff(double x)
{
    return x*x;
}
int k;
int main()
{
    scanf("%lf%lf",&x[0],&c);
    while (1)
    {
        w[k]=1/(1-ff(x[k]));
        x[k+1]=(1-w[k])*x[k]+w[k]*f(x[k]);
        if (abs(x[k+1]-x[k])<c ) 
        {
            printf("%.19lf\n",x[k+1]);
            break;
        }else k=k+1;
    }
}
View Code

2.10-牛顿法求根

 1 #include<bits./stdc++.h>
 2 using namespace std;
 3 double d,c,A,B,x[1000],w[1000];
 4 double f(double x)
 5 {
 6     return x*x*x-3.0*x+1;
 7 } 
 8 double ff(double x)
 9 {
10     return 3*x*x-3;
11 }
12 int k;
13 int main()
14 {
15     scanf("%lf%lf",&x[0],&c);
16     k=0;
17     while(1)
18     {
19         x[k+1]=x[k]-f(x[k])/ff(x[k]);
20         if (abs(x[k+1] -x[k])  < c) 
21         {
22             printf("%.19lf\n",x[k+1]);
23             break;
24         }else k=k+1;
25     }
26 }
View Code
posted @ 2019-06-13 22:30  口香糖万岁  阅读(1545)  评论(0)  编辑  收藏  举报