实习题11用两种不容的顺序计算644834.1100012≈∑=-n n，分析误差的变化（1）顺序计算源代码：#include<stdio.h>#include<math.h>void main(){double sum=0;int n=1;while(1){sum=sum+(1/pow(n,2));if(n%1000==0)printf("sun[%d]=%-30f",n,sum);if(n>=10000)break;n++;}printf("sum[%d]=%f\n",n,sum);}结果：（2）逆序计算源代码：#include<stdio.h>#include<math.h>void main(){double sum=0;int n=10000;while(1){sum=sum+(1/pow(n,2));if(n%1000==0)printf("sum[%d]=%-30f",n,sum);if(n<=1)break;n--;}printf("sum[%d]=%f\n",n,sum);}结果：2已知连分数))//(.../(32211nnbaabababf++++=利用下面的方法计算f:11)0,...,2,1(,dfnnidabdbdiiiinn=--=+==++写一个程序，读入n,nnba,,计算并打印f源代码：#include<stdio.h>#include<math.h>void main(){int i=0,n;float a[1024],b[1024],d[1024];printf("please input n,n=");scanf("%d",&n);printf("\nplease input a[1] to a[n]:\n");for(i=1;i<=n;i++){printf("a[%d]=",i);scanf("%f",&a[i]);}printf("\nplease input b[0] to b[n]:\n");for(i=0;i<=n;i++){printf("b[%d]=",i);scanf("%f",&b[i]);}d[n]=b[n];for(i=n-1;i>=0;i--)d[i]=b[i]+a[i+1]/d[i+1];printf("\nf=%f\n",d[0]);}结果：3给出一个有效的算法和一个无效的算法计算积分⎰=+=10)10,...1,0(14n dx x x y nn 源代码：#include<stdio.h>#include<math.h>main(){double y_0=(1/4.0)*log(5),y_1;double y_2=(1.0/55.0+1.0/11.0)/2,y_3;int n=1,m=10;printf("有效算法输出结果：\n");printf("y[0]=%-20f",y_0);while(1){y_1=1.0/(4*n)+y_0/(-4.0);printf("y[%d]=%-20f",n,y_1);if(n>=10)break;y_0=y_1;n++;if(n%3==0)printf("\n");}printf("\n 无效算法的输出结果：\n");printf("y[10]=%-20f",y_2);while(1){y_3=1.0/n-4.0*y_2;printf("y[%d]=%-20f",m-1,y_3);if(m<=1)break;y_2=y_3;m--;if(m%2==0) printf("\n");}}结果：4设∑=-=N j N j S 2211，已知其精确值为)11123(21+--N N （1）编制按从小到大顺序计算N S 的程序（2）编制按从小达到的顺序计算N S 的程序（3）按两种顺序分别计算30000100001000,,S S S ,并指出有效位数源代码：#include<stdio.h>main(){int N;double SN[30000];SN[30000]=(3.0/2.0-1.0/30000.0-1/30001.0)/2.0;for(N=30000;N>=2;N--)SN[N-1]=SN[N]-1.0/(N*N-1);printf("从大到小顺序计算： \nSN[1000]=%f\nSN[10000]=%f\nSN[30000]=%f\n",SN[1000],SN[10000],SN[30000]); SN[2]=(3.0/2-1.0/2.0-1/3.0)/2.0;for(N=3;N<=30000;N++)SN[N]=SN[N-1]+1.0/(N*N-1);printf("从小到大顺序计算： \nSN[1000]=%f\nSN[10000]=%f\nSN[30000]=%f\n",SN[1000],SN[10000],SN[30000]); }结果：实习题21.用牛顿法求下列方程的根02=-x e x01=-x xe02lg =-+x x源代码：#include <stdio.h>#include <math.h>typedef float (*p)(float );float ff1(float x){return x*x-exp(x);}float ff2(float x){return x*exp(x)-1;}float ff3(float x){return log(x)+x-2;}float answer(float(*p)(float)){int k=2;float m=1,n=-1,x2,a,b,c;if (p==ff3)n=2;printf("x[0] = %.4f, x[1] = %.4f, ",m,n);while (1){if (fabs(m-n)<1e-4) break;a=p(n)*(n-m);b=p(n)-p(m);c=a/b;x2=n-c;m = n;n = x2;printf("x[%d] = %.4f, ",k,x2);k++;if (k%3==0) printf("\n");}if (k%3!=0) printf("\n");printf("iteration times: %d, roots: %.4f\n ",k-2,n);return 0;}main(){printf("x*x-exp(x),\n");answer(ff1);printf("x*exp(x)-1,\n");answer(ff2);printf("lg(x)+x-2,\n");answer(ff3);return 0;}结果：2.编写一个割线法的程序，求解上述各方程源代码：#include<stdio.h>#include<math.h>float gexian(float,float);float f(float);main(){int i,j;float x1=2.2;float x2=2,x3;scanf("%d",&i);if(i==1)printf("%f",x1);else if(i==2)printf("%f",x2);else{for(j=3;j<=i;j++){x3=gexian(x1,x2);x1=x2;x2=x3;}printf("%f",gexian(x1,x2));}}float f(float x){return (x*x-exp(x));}float gexian(float x1,float x2){return (x2-(f(x2)/(f(x2)-f(x1)))*(x2-x1));}结果：实习题3 1用列主元消去法解下列方程组；⎪⎪⎩⎪⎪⎨⎧=++=-++--=+--=--+43443233312)1(421432143214321xxxxxxxxxxxxxxx⎪⎪⎩⎪⎪⎨⎧=++--=++-=-+--=-+-4341220332282)2(432132143214321xxxxxxxxxxxxxxx源程序：#include<stdio.h>#include<math.h>void ColPivot(float*,int,float[]);void ColPivot(float*c,int n,float x[]){int i,j,t,k;float p;for(i=0;i<=n-2;i++){k=i;for(j=i+1;j<=n-1;j++)if(fabs(*(c+j*(n+1)+i))>(fabs(*(c+k*(n+1)+i))))k=j;if(k!=i)for(j=i;j<=n;j++){p=*(c+i*(n+1)+j);*(c+i*(n+1)+j)=*(c+k*(n+1)+j);*(c+k*(n+1)+j)=p;}for(j=i+1;j<=n-1;j++){p=(*(c+j*(n+1)+i))/(*(c+i*(n+1)+i));for(t=i;t<=n;t++)*(c+j*(n+1)+t)-=p*(*(c+i*(n+1)+t));}}for(i=n-1;i>=0;i--){for(j=n-1;j>=i+1;j--)(*(c+i*(n+1)+n))-=x[j]*(*(c+i*(n+1)+j));x[i]=*(c+i*(n+1)+n)/(*(c+i*(n+1)+i));}}void main(){int i;float x[4];float c[4][5]={1,1,0,3,4,2,1,-1,1,1,3,-1,-1,3,-3,-1,2,3,-1,4};ColPivot(c[0],4,x);for(i=0;i<=3;i++)printf("x[%d]=%f\n",i,x[i]);}结果：第（1）题第（2）题2、源代码：#include<stdio.h>void main(){float x[4];int i;float a[4][5]={48,-24,0,-12,4,-24,24,12,12,4,0,6,20,2,-2,-6,6,2,16,-2};void DirectLU(float*,int,float[]);DirectLU(a[0],4,x);for(i=0;i<=3;i++)printf("x[%d]=%f\n",i,x[i]);}void DirectLU(float*u,int n,float x[]){int i,r,k;for(r=0;r<=n-1;r++){for(i=r;r<=n;i++)for(k=0;k<=r-1;k++)*(u+r*(n+1)+i)-=*(u+r*(n+1)+k)*(*(u+k*(n+1)+i));for(i=r+1;i<=n-1;i++){for(k=0;k<=r-1;k++)*(u+i*(n+1)+r)-=*(u+i*(n+1)+k)*(*(u+k*(n+1)+r));*(u+i*(n+1)+r)/=*(u+r*(n+1)+r);}}for(i=n-1;i>=0;i--){for(r=n-1;r>=i+1;r--)*(u+i*(n+1)+n)-=*(u+i*(n+1)+r)*x[r];x[i]=*(u+i*(n+1)+n)/(*(u+i*(n+1)+i));}}实习题41、源代码：#include<stdio.h>float Lagrange(float x[],float y[],float xx,int n) //n为（n+1）次插值；{int i,j;float *a,yy=0;a=new float[n];for(i=0;i<=n-1;i++){a[i]=y[i];for(j=0;j<=n-1;j++)if(j!=i)a[i]*=(xx-x[j])/(x[i]-x[j]);yy+=a[i];}delete a;return yy;}void main(){float x[5]={-3.0,-1.0,1.0,2.0,3.0};float y[5]={1.0,1.5,2.0,2.0,1.0};float xx1=-2,xx2=0,xx3=2.75,yy1,yy2,yy3;yy1=Lagrange(x,y,xx1,3);yy2=Lagrange(x,y,xx2,3);yy3=Lagrange(x,y,xx3,3);printf("x1=%-20f,y1=%f\n",xx1,yy1);printf("x2=%-20f,y2=%f\n",xx2,yy2);printf("x3=%-20f,y3=%f\n",xx3,yy3);}结果：2、源代码：#include<stdio.h>float Lagrange(float x[],float y[],float xx,int n) //n为（n+1）次插值；{int i,j;float *a,yy=0;a=new float[n];for(i=0;i<=n-1;i++){a[i]=y[i];for(j=0;j<=n-1;j++)if(j!=i)a[i]*=(xx-x[j])/(x[i]-x[j]);yy+=a[i];}delete a;return yy;}void main(){float x[6]={0.30,0.42,0.50,0.58,0.66,0.72};float y[6]={1.04403,1.08462,1.11803,1.15603,1.19817,1.23223};float xx1=0.46,xx2=0.55,xx3=0.60,yy1,yy2,yy3;yy1=Lagrange(x,y,xx1,6);yy2=Lagrange(x,y,xx2,6);yy3=Lagrange(x,y,xx3,6);printf("x1=%-20f,y1=%f\n",xx1,yy1);printf("x2=%-20f,y2=%f\n",xx2,yy2);printf("x3=%-20f,y3=%f\n",xx3,yy3);}结果：源代码：#include<stdio.h>#define N 3void Difference(float y[],float f[4][4],int n){int k,i;f[0][0]=y[0];f[1][0]=y[1];f[2][0]=y[2];f[3][0]=y[3];for(k=1;k<=n;k++)for(i=0;i<=(N-k);i++)f[i][k]=f[i+1][k-1]-f[i][k-1];return;}void main(){int i,k=1;float a,b=1,m=21.4,t=1.4,f[4][4]={0};float x[5]={20,21,22,23,24};float y[5]={1.30103,1.32222,1.34242,1.36173,1.38021};Difference(y,f,N);a=f[0][0];for(i=1;i<=N;i++){k=k*i;b=b*(t-i+1);a=a+b*f[0][i]/k;}printf("x(k)\n");for (i=0;i<=4;i++)printf( "%-20f",x[i]);printf("\ny(k)\n");for (i=0;i<=4;i++)printf("%-20f",y[i]);for(k=1;k<=3;k++){printf("\nF(%d)\n ",k);for(i=0;i<=(3-k);i++){printf("%-20f",f[i][k]);}}printf ("\n");printf("f(%f)=%-20f",m,a);printf ("\n");结果：实习题52、源代码：#include<stdio.h>#include<math.h>void main(){int i,n;float a[2];float x[15]={1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8},z[15];floaty[15]={33.4,79.50,122.65,159.05,189.15,214.15,238.65,252.50,267.55,280.50,296.65,30 1.40,310.40,318.15,325.15};for(n=0;n<=14;n++) //增加了数组z；{z[n]=log(y[n]/x[n]);}void Approx(float[],float[],int,int,float[]);Approx(x,z,15,1,a); //变成一次拟合；//for(i=0;i<=1;i++)//printf("a[%d]=%f\n",i,a[i]);printf("a=exp(a[0])=%f\n",exp(a[0]));printf("b=-a[1]=%f\n",-a[1]); }void Approx(float x[],float y[],int m,int n,float a[]){int i,j,t;float *c=new float[(n+1)*(n+2)];float power(int,float);void ColPivot(float *,int,float[]);for(i=0;i<=n;i++){for(j=0;j<=n;j++){*(c+i*(n+2)+j)=0;for(t=0;t<=m-1;t++)*(c+i*(n+2)+j)+=power(i+j,x[t]);}*(c+i*(n+2)+n+1)=0;for(j=0;j<=m-1;j++)*(c+i*(n+2)+n+1)+=y[j]*power(i,x[j]);}ColPivot(c,n+1,a);delete c;}void ColPivot(float *c,int n,float x[]){int i,j,t,k;float p;for(i=0;i<=n-2;i++){k=i;for(j=i+1;j<=n-1;j++)if(fabs(*(c+j*(n+1)+i))>(fabs(*(c+k*(n+1)+i))))k=j;if(k!=i)for(j=i;j<=n;j++){p=*(c+i*(n+1)+j);*(c+i*(n+1)+j)=*(c+k*(n+1)+j);*(c+k*(n+1)+j)=p;}for(j=i+1;j<=n-1;j++){p=(*(c+j*(n+1)+i))/(*(c+i*(n+1)+i));for(t=i;t<=n;t++)*(c+j*(n+1)+t)-=p*(*(c+i*(n+1)+t));}}for(i=n-1;i>=0;i--){for(j=n-1;j>=i+1;j--)(*(c+i*(n+1)+n))-=x[j]*(*(c+i*(n+1)+j));x[i]=*(c+i*(n+1)+n)/(*(c+i*(n+1)+i));}float power(int i,float v){float a=1;while(i--)a*=v;return a;}结果：实习题61、源代码：（1）#include<stdio.h>#include<math.h>float Cotes(float(*f)(float),float a,float b,int n){int k;float c,c1=0,c2,c3,c4;float h=(b-a)/n;c2=(*f)(a+h/4);c3=(*f)(a+h/2);c4=(*f)(a+3*h/4);for(k=1;k<=n-1;k++){c1+=(*f)(a+k*h);c2+=(*f)(a+k*h+h/4);c3+=(*f)(a+k*h+h/2);c4+=(*f)(a+k*h+3*h/4);}c=h/90*(7*((*f)(a)+(*f)(b))+14*c1+32*c2+12*c3+32*c4);return c;}float f(float x){return 1/sqrt(1+x*x*x);}void main()int i,n=4;float c;for(i=0;i<=4;i++){c=Cotes(f,0,1,n);printf("C(%d)=%f\n",n,c);n*=2;}}（2）#include<stdio.h>#include<math.h>float Cotes(float(*f)(float),float a,float b,int n){int k;float c,c1=0,c2,c3,c4;float h=(b-a)/n;c2=(*f)(a+h/4);c3=(*f)(a+h/2);c4=(*f)(a+3*h/4);for(k=1;k<=n-1;k++){c1+=(*f)(a+k*h);c2+=(*f)(a+k*h+h/4);c3+=(*f)(a+k*h+h/2);c4+=(*f)(a+k*h+3*h/4);}c=h/90*(7*((*f)(a)+(*f)(b))+14*c1+32*c2+12*c3+32*c4);return c;}float f(float x){ // return 1/sqrt(1+x*x*x);if (x==0)return 1;else return sin(x)/x;}void main(){int i,n=4;float c;for(i=0;i<=4;i++){// c=Cotes(f,0,1,n);c=Cotes(f,0,5,n);printf("C(%d)=%f\n",n,c);n*=2;}}结果：（1）（2）实习题7 一、改进欧拉法1、#include<stdio.h>#include<iostream>double Adams (double (*f)(double x, double y),double x0,double y0,double h,int N) {for(int n=0; n<N; ++n) {double x1=x0+h;double yp=y0+h*f(x0,y0);double y1=y0+h*f(x1,yp);y1=(yp+y1)/2.0;printf("ty=%f\n",y1);x0=x1;y0=y1;}}int main(void){double f(double x, double y); double x0=0,y0=0;double x,y,step;long i;step=0.1;Adams(f,x0,y0,0.1,10);}double f(double x, double y) {double r;r=x*x+y*y;return(r);}2、int main(void){double f(double x, double y); double x0=0,y0=0;double x,y,step;long i;step=0.1;Adams(f,x0,y0,0.1,10);}double f(double x, double y) {double r;r=1/(1+y*y);return(r);}3、int main(void){double f(double x, double y); double x0=0,y0=1;double x,y,step;long i;step=0.1;Adams(f,x0,y0,0.1,10);}double f(double x, double y) {double r;r=y-2*x/y;return(r);}四阶龙格库塔法1、#include<iostream>#include<stdio.h>using namespace std;//f为函数的入口地址，x0、y0为初值，xn为所求点，step为计算次数double Runge_Kuta( double (*f)(double x, double y), double x0, double y0, double xn, long step ){double k1,k2,k3,k4,result;double h=(xn-x0)/step;if(step<=0)return(y0);if(step==1){k1=f(x0,y0);k2=f(x0+h/2, y0+h*k1/2);k3=f(x0+h/2, y0+h*k2/2);k4=f(x0+h, y0+h*k3);result=y0+h*(k1+2*k2+2*k3+k4)/6;}else{double x1,y1;x1=xn-h;y1=Runge_Kuta(f, x0, y0, xn-h,step-1);k1=f(x1,y1);k2=f(x1+h/2, y1+h*k1/2);k3=f(x1+h/2, y1+h*k2/2);k4=f(x1+h, y1+h*k3);result=y1+h*(k1+2*k2+2*k3+k4)/6;}printf("%lg\n",result);return(result);}int main(void){double f(double x, double y);double x0=0,y0=0,k;double x,y,step;long i;step=0.1;cout.precision(10);//for(i=0;i<=10;i++)//{// x=x0+i*step;// cout<<x<<" - "<<Runge_Kuta(f,x0,y0,x,i)<<endl; //}//cout<<Runge_Kuta(f,x0,y0,1,10);k= Runge_Kuta(f,x0,y0,1,10);}double f(double x, double y){double r;r=x*x+y*y;return(r);}2、int main(void){double f(double x, double y);double x0=0,y0=0,k;double x,y,step;long i;step=0.1;cout.precision(10);//for(i=0;i<=10;i++)//{// x=x0+i*step;// cout<<x<<" - "<<Runge_Kuta(f,x0,y0,x,i)<<endl; //}//cout<<Runge_Kuta(f,x0,y0,1,10);k= Runge_Kuta(f,x0,y0,1,10);}double f(double x, double y){double r;r=1/（1+y*y）;return(r);}3、int main(void){double f(double x, double y);double x0=0,y0=0,k;double x,y,step;long i;step=0.1;cout.precision(10);//for(i=0;i<=10;i++)//{// x=x0+i*step;// cout<<x<<" - "<<Runge_Kuta(f,x0,y0,x,i)<<endl; //}//cout<<Runge_Kuta(f,x0,y0,1,10);k= Runge_Kuta(f,x0,y0,1,10);}double f(double x, double y){double r;r=1/（1+y*y）;return(r);}二、int main(void){double f(double x, double y);double x0=0,y0=1,k;double x,y,step;long i;step=0.1;cout.precision(10);//for(i=0;i<=10;i++)//{// x=x0+i*step;// cout<<x<<" - "<<Runge_Kuta(f,x0,y0,x,i)<<endl;//}//cout<<Runge_Kuta(f,x0,y0,1,10);k= Runge_Kuta(f,x0,y0,1,10);}double f(double x, double y){double r;r=y-2*x/y;return(r);}三、1、void Runge_Kutta(float(*f)(float x,float y),float a,float b,float y0,int N,float yy[]) {float x=a,y=y0,K1,K2,K3,K4;float h=(b-a)/N;int i;for(i=1;i<=3;i++){K1=(*f)(x,y);K2=(*f)(x+h/2,y+h*K1/2);K3=(*f)(x+h/2,y+h*K2/2);K4=(*f)(x+h/2,y+h*K3);y=y+h*(K1+2*K2+2*K2+2*K3+K4)/6;x=a+i*h;yy[i-1]=y;}}void Adams (float a,float b,int N,float(*f)(float x,float y),float y0){int i;float y1,y2,y,yp,yc,yy[3],h,x;printf("x[0]=%f\ty[0]=%f\n",a,y0);Runge_Kutta(f,a,b,y0,N,yy);y1=yy[0];y2=yy[1];y=yy[2];h=(b-a)/N;for(i=1;i<=3;i++)printf("x[%d]=%f\ty[%d]=%f\n",i,a+i*h,i,yy[i-1]);for(i=3;i<N;i++){x=a+i*h;yp=y+h*(55*f(x,y)-59*f(x-h,y2)+37*f(x-2*h,y1)-9*f(x-3*h,y0))/24;yc=y+h*(9*(*f)(x+h,yp)+19*(*f)(x,y)-5*(*f)(x-h,y2)+(*f)(x-2*h,y1))/24;printf("x[%d]=%f\ty[%d]=%f\n",i+1,x+h,i+1,yc);y0=y1;y1=y2;y2=y;y=yc;}}float f(float x,float y){return y*y;void Runge_Kutta(float(*f)(float x,float y),float a,float,float b,float y0,int N,float yy[]); void Adams (float a,float b,int N,float(*f)(float x,float y),float y0);float f(float x,float y);int main (void){float a=0,b=1,y0=1;int N=10;Adams(a,b,N,f,y0);}2、#include<stdio.h>void Runge_Kutta(float(*f)(float x,float y),float a,float b,float y0,int N,float yy[]){float x=a,y=y0,K1,K2,K3,K4;float h=(b-a)/N;int i;for(i=1;i<=3;i++){K1=(*f)(x,y);K2=(*f)(x+h/2,y+h*K1/2);K3=(*f)(x+h/2,y+h*K2/2);K4=(*f)(x+h/2,y+h*K3);y=y+h*(K1+2*K2+2*K2+2*K3+K4)/6;x=a+i*h;yy[i-1]=y;}}void Adams (float a,float b,int N,float(*f)(float x,float y),float y0)int i;float y1,y2,y,yp,yc,yy[3],h,x;printf("x[0]=%f\ty[0]=%f\n",a,y0);Runge_Kutta(f,a,b,y0,N,yy);y1=yy[0];y2=yy[1];y=yy[2];h=(b-a)/N;for(i=1;i<=3;i++)printf("x[%d]=%f\ty[%d]=%f\n",i,a+i*h,i,yy[i-1]);for(i=3;i<N;i++){x=a+i*h;yp=y+h*(55*f(x,y)-59*f(x-h,y2)+37*f(x-2*h,y1)-9*f(x-3*h,y0))/24;yc=y+h*(9*(*f)(x+h,yp)+19*(*f)(x,y)-5*(*f)(x-h,y2)+(*f)(x-2*h,y1))/24;printf("x[%d]=%f\ty[%d]=%f\n",i+1,x+h,i+1,yc);y0=y1;y1=y2;y2=y;y=yc;}}float f(float x,float y){return 0.1*(x*x*x+y*y);}void Runge_Kutta(float(*f)(float x,float y),float a,float,float b,float y0,int N,float yy[]); void Adams (float a,float b,int N,float(*f)(float x,float y),float y0);float f(float x,float y);int main (void){float a=0,b=1,y0=1;int N=10;Adams(a,b,N,f,y0);}。