1.1-有效数字丢失现象观察
1 #include2 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 }
1.2-n次插值的Lagrange形式
1 #include2 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 }
1.3-n次插值Newton形式
1 #include2 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 }
1.4-三次样条插值
1 #include2 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 }
1.5-单变量数据拟合（最小二乘法）
1 #include2 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 }
1.6-自动选取步长梯形法
1 #include2 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 }
1.7-Romberg 求积法
1 #include2 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])
2.1-顺序高斯消去法
1 #include2 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<
2.2-列主元高斯消去法
1 #include2 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<
2.3-全主元高斯消去法
1 #include2 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<
2.4-LU 直接分解法
1 #include2 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<
2.5-Jacobi 迭代算法
1 #include2 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])
2.6-Seidel 迭代算法
1 #include2 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])
2.7-松弛迭代算法
1 #include2 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])
2.8-对分法求根
1 #include2 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)
2.9-松弛法迭代求根
#includeusing 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])
2.10-牛顿法求根
1 #include2 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 }