программы_вычмат. Контрольная работа 1 по дисциплине Вычислительная математика учебное пособие Мицель А. А. Вычислительная математика
 Скачать 451.5 Kb.
|
1 2 Листинг вычислительных алгоритмов Листинг программы задания 1 //--------------------------------------------------------------------------- #include #include #include #pragma hdrstop //--------------------------------------------------------------------------- //функция, для которой ищем корни f(x)=0 double f(float x) { return 2*x*x-5-pow(2,x); } //поиск корней методом золотого сечения double zolotoe(double a, double b, double eps, int &n, double &alpha, double &toch) { double c,d,x0,x1,x2; //точки на отрезке double lambda=(sqrt(5)+1)/2.0; n=0; //число итераций do { d=a+(b-a)/lambda; c=b-d+a; if (f(a)*f(d)<0) //определяем, на какой фрагмент отрезка попадает корень b=d; else a=c; n++; x0=x1; x1=x2; x2=(a+b)/2; //середина отрезка } while (((b-a)>=2*eps)||fabs(f(x2))>=eps); alpha=fabs((x2-x1)/(x1-x0)); //параметр сходимости toch=(b-a)/2; //точность корня return x2; } int main(int argc, char* argv[]) { double a, b, eps; printf("a="); scanf("%lf", &a); printf("b="); scanf("%lf", &b); printf("eps="); scanf("%lf", &eps); double alpha, toch; int n; double x=zolotoe(a,b,eps,n,alpha,toch); printf("x=%f \n\r",x,toch,eps); printf("tochnost: %e (eps=%e)\n\r", toch, eps); printf("f(x)=%f \n\r",f(x)); printf("Count of iterations N= %d \n\r",n); printf("Count of evaluating f(x) %d \n\r",n*3); printf("alpha=%e", alpha); getch(); return 0; } Листинг программы задания 2.1 //--------------------------------------------------------------------------- #pragma hdrstop #include #include #include #include #include FILE *FileOutput; //--------------------------------------------------------------------------- #pragma argsused //меняет строки матрицы местами void swap(double *a, double *b, int n) { double temp; for (int i=0;i temp=a[i]; a[i]=b[i]; b[i]=temp; } } //вывод матрицы в файл void PrintMatrix(double **a,int n,char* s, char* c = "%lf ") { fprintf(FileOutput,"\n\n%s:\n",s); for(int i=0;i fprintf(FileOutput,"\n"); for(int j=0;j fprintf(FileOutput,c,a[i][j]); } } } //вывод вектора в файл void PrintVector(double *a,int n,char* s, char *c = "%lf ") { fprintf(FileOutput,"\n\n%s:\n",s); for(int i=0;i fprintf(FileOutput,c,a[i]); } } //норма вектора double normal_vector(double *epsilon,int n) { double normal=0; for(int i=0;i normal+=epsilon[i]*epsilon[i]; } normal=sqrt(normal); return normal; } //норма матрицы double normal_matrix(double **e,int n) { double normal=0; for(int i=0;i for(int j=0;j normal+=e[i][j]*e[i][j]; } } normal=sqrt(normal); return normal; } //вычитание векторов double * MinusVector(double *a,double *b,int n) { double *c; c=(double*)malloc(n*sizeof(double)); for(int i=0;i c[i]=b[i]-a[i]; } return c; } //решение СЛАУ методом Гаусса double* Gauss(double **a,int n) { double *x; x=(double*)malloc(n*sizeof(double)); //прямой ход метода for (int i=0; i int k = i; for (int j=i+1; j if (fabs (a[j][i]) > fabs (a[k][i])) k = j; swap (a[i],a[k],n); for (int j=i+1; j a[i][j] /= a[i][i]; a[i][i]=1; for (int j=i+1; j double koef=a[j][i]; for (int k=0; k a[j][k] -= a[i][k] * koef; } } //обратный ход for (int i=n-1; i>=0; i--) { for (int j=i+1; j a[i][j] /= a[i][i]; for (int j=0; j double koef=a[j][i]; for (int k=0; k a[j][k] -= a[i][k] * koef; } } //выделяем обратную матрицу for (int i=0;i x[i]=a[i][n]; return x; } //вычисление невязки double* Verification_SLAU(double **a,double *x, int n) { double *e; e=(double*)malloc(n*sizeof(double)); for(int i=0;i e[i]=0; for(int j=0;j e[i]=a[i][j]*x[j]+e[i]; } e[i]=a[i][n]-e[i]; } return e; } int main(int argc, char* argv[]) { FILE *FileInput; int n; double **a,*x; double eps,normal; FileInput=fopen("input.txt","r"); fscanf(FileInput, "%lf", &eps); fscanf(FileInput, "%d", &n); a=(double**)malloc(n*sizeof(double*)); for (int i=0;i a[i]=(double*)malloc((n+1)*sizeof(double)); for(int j=0;j<=n;j++) { fscanf(FileInput,"%lf",&a[i][j]); } } fclose(FileInput); FileOutput=fopen("output.txt","w"); double **X; x=Gauss(a,n); PrintVector(x,n,"Decision SLAU"); x=Verification_SLAU(a,x,n); PrintVector(x,n,"Verification SLAU","%e "); double d=normal_vector(x,n); fprintf(FileOutput,"\n\nNorma\n\n%e",d); free(x); return 0; } Листинг программы задания 2.2 #pragma hdrstop #include #include #include #include #include FILE *FileOutput; //--------------------------------------------------------------------------- #pragma argsused //меняет строки матрицы местами (точнее, любые массивы размера n) void swap(double *a, double *b, int n) { double temp; for (int i=0;i temp=a[i]; a[i]=b[i]; b[i]=temp; } } //вычисление определителя методом Гаусса double Gauss(double **a,int n,double eps) { double det = 1; for (int i=0; i int k = i; for (int j=i+1; j if (fabs (a[j][i]) > fabs (a[k][i])) k = j; //если весь столбец нулевой, //то определитель также будет равен 0 if (fabs(a[k][i]) < eps) { det = 0; break; } swap (a[i], a[k],n); if (i != k) det = -det; det *= a[i][i]; for (int j=i+1; j a[i][j] /= a[i][i]; for (int j=0; j if (j != i && fabs (a[j][i]) > eps) for (int k=i+1; k a[j][k] -= a[i][k] * a[j][i]; } return det; } int main(int argc, char* argv[]) { FILE *FileInput; int n; double **a; double eps; FileInput=fopen("input.txt","r"); fscanf(FileInput, "%lf", &eps); fscanf(FileInput, "%d", &n); a=(double**)malloc(n*sizeof(double*)); for (int i=0;i a[i]=(double*)malloc(n*sizeof(double)); for(int j=0;j fscanf(FileInput,"%lf",&a[i][j]); } } fclose(FileInput); FileOutput=fopen("output.txt","w"); double d=Gauss(a,n,eps); fprintf(FileOutput,"\n\det=%f",d); free(a); } Листинг программы задания 2.3 #pragma hdrstop #include #include #include #include #include FILE *FileOutput; //--------------------------------------------------------------------------- #pragma argsused //меняет строки матрицы местами (точнее, любые массивы размера n) void swap(double *a, double *b, int n) { double temp; for (int i=0;i temp=a[i]; a[i]=b[i]; b[i]=temp; } } //вывод матрицы в файл void PrintMatrix(double **a,int n,char* s, char* c = "%f ") { fprintf(FileOutput,"\n\n%s:\n",s); for(int i=0;i fprintf(FileOutput,"\n"); for(int j=0;j fprintf(FileOutput,c,a[i][j]); } } } //обратная матрица методом Гаусса double** Gauss(double** a_gl,int n,double eps) { double **x, **a; //чтобы не портить матрицу //копируем ее в новую а x=(double**)malloc(n*sizeof(double)); a=(double**)malloc(2*n*sizeof(double)); for (int i=0;i x[i]=(double*)malloc(n*sizeof(double)); a[i]=(double*)malloc(2*n*sizeof(double)); a[i][n+i]=1; for (int j=0;j a[i][j]=a_gl[i][j]; } //прямой ход метода for (int i=0; i int k = i; for (int j=i+1; j if (fabs (a[j][i]) > fabs (a[k][i])) k = j; swap (a[i], a[k],2*n); for (int j=i+1; j<2*n; j++) a[i][j] /= a[i][i]; a[i][i]=1; for (int j=i+1; j double koef=a[j][i]; for (int k=0; k<2*n; k++) a[j][k] -= a[i][k] * koef; } } //обратный ход for (int i=n-1; i>=0; i--) { for (int j=i+1; j<2*n; j++) a[i][j] /= a[i][i]; for (int j=0; j double koef=a[j][i]; for (int k=0; k<2*n; k++) a[j][k] -= a[i][k] * koef; } } //выделяем обратную матрицу for (int i=0;i for (int j=0;j x[i][j]=a[i][j+n]; return x; } //умножение двух квадратных матриц одной размерности double** MultMatrix(double **a, double **b, int n) { double **y; y=(double**)malloc(n*sizeof(double)); for (int i=0;i y[i]=(double*)malloc(n*sizeof(double)); for (int i=0;i for (int j=0;j y[i][j]=0; for (int k=0;k y[i][j]+=a[i][k]*b[k][j]; } return y; } int main(int argc, char* argv[]) { FILE *FileInput; int n; double **a,**x; double eps; FileInput=fopen("input.txt","r"); fscanf(FileInput, "%lf", &eps); fscanf(FileInput, "%d", &n); a=(double**)malloc(n*sizeof(double*)); for (int i=0;i a[i]=(double*)malloc((2*n)*sizeof(double)); for(int j=0;j fscanf(FileInput,"%lf",&a[i][j]); } } fclose(FileInput); FileOutput=fopen("output.txt","w"); x=Gauss(a,n,eps); PrintMatrix(x,n,"A^-1="); x=MultMatrix(a,x,n); for (int i=0;i x[i][i]-=1; } PrintMatrix(x,n,"A*A^-1 - E=","%e "); return 0; } //--------------------------------------------------------------------------- Листинг программы задания 3 //--------------------------------------------------------------------------- #pragma hdrstop #include #include #include #include //--------------------------------------------------------------------------- #pragma argsused //функция, которую интерполируем double y(double x) { return 1-2*tan(x); } // интерполирование Лагранжа double Lagrang(double x, int n, double h) { int b=x/h; int i_st, i_fin; if (n>11) n=11; //больше по условию не задано if ((b-n/2)>0) { i_st=b-n/2; i_fin=b+(n-n/2); if (i_fin>11) { i_fin=11; i_st=11-n; } } else { i_st=0; i_fin=n; } //вычисление полинома double s=0; for (int i=i_st;i double p=1; for (int j=i_st;j if (i!=j) p*=(x-j*h)/(i*h-j*h); s+=y(i*h)*p; } return s; } int main(int argc, char* argv[]) { int n; double h=0.02*M_PI; double hnew=h/2; //новый шаг - половина от старого printf("Polynomial degree="); scanf("%i", &n); printf(" x y Polynom delta\n\r"); for (double x=0; x<0.2*M_PI+1E-7; x+=hnew) { double L=Lagrang(x,n,h); printf("%f %f %f %e\n\r", x,y(x),L,fabs(L-y(x))); } getch(); return 0; } Листинг программы задания 4 //--------------------------------------------------------------------------- #pragma hdrstop #include #include #include #include //--------------------------------------------------------------------------- #pragma argsused //функция, которую интерполируем double y(double x) { return 3.0/tan(x); } //полином Лагранжа double Lagrang(double x, int n, int max, double h) { if (n>max) n=max; //больше точек по условию не задано int b=(x-1)/h; int i_st, i_fin; if ((b-n/2)>0) { i_st=b-n/2; i_fin=b+(n-n/2); if (i_fin>max) { i_fin=max; i_st=max-n; } } else { i_st=0; i_fin=n; } //вычисление полинома double s=0; for (int i=i_st;i double p=1; for (int j=i_st;j if (i!=j) p*=(x-(0.3*M_PI+j*h))/((0.3*M_PI+i*h)-(0.3*M_PI+j*h)); s+=y((0.3*M_PI+i*h))*p; } return s; } //первая производная double diff(double x, int n, int max, double h) { int b=(x-0.3*M_PI)/h; if (b==0) return (-3*Lagrang(x,n,max,h)+4*Lagrang(x+h,n,max,h)-Lagrang(x+2*h,n,max,h))/(2*h); if (b==max-1) return (Lagrang(x-2*h,n,max,h)-4*Lagrang(x-h,n,max,h)+3*Lagrang(x,n,max,h))/(2*h); return (Lagrang(x+h,n,max,h)-Lagrang(x-h,n,max,h))/(2*h); } //вторая производная double diff2(double x, int n, int max, double h) { int b=(x-0.3*M_PI)/h; if (b==0) return (2*Lagrang(x,n,max,h)-5*Lagrang(x+h,n,max,h)+4*Lagrang(x+2*h,n,max,h)-Lagrang(x+3*h,n,max,h))/(h*h); if (b==max-1) return (-Lagrang(x-3*h,n,max,h)+4*Lagrang(x-2*h,n,max,h)-5*Lagrang(x-h,n,max,h)+2*Lagrang(x,n,max,h))/(h*h); return (Lagrang(x-h,n,max,h)-2*Lagrang(x,n,max,h)+Lagrang(x+h,n,max,h))/(h*h); } int main(int argc, char* argv[]) { double h=0.02*M_PI; double hnew=h/2; int n; printf("Polinomial degree="); scanf("%i", &n); printf(" x diff diff err diff2 diff2 err \n\r"); for (double x=0.3*M_PI; x<0.5*M_PI+1E-7; x+=hnew) { printf("%f %f %e %f %e\n\r", x,diff(x,n,11,h),fabs(diff(x,n,11,h)+3/pow(sin(x),2)),diff2(x,n,11,h),fabs(diff2(x,n,11,h)-6*cos(x)/pow(sin(x),3))); } getch(); return 0; } Листинг программы задания 5 //--------------------------------------------------------------------------- #pragma hdrstop #include #include #include #include //--------------------------------------------------------------------------- //функция, которую нужно интегрировать double f(double x) { return 1/(x*x+4); } //метод Симпсона double Simpson(double a, double b, double eps, int &n, double &er) { n=2; double I0,I1=(f(a)+4*f((a+b)/2)+f(b))*(b-a)/3; do { n*=2; double sum=f(a)+f(b); double h=1.0*(b-a)/n; double x=a; for (int i=1;i x+=h; if (i%2==1) sum+=4*f(x); else sum+=2*f(x); } I0=I1; I1=sum*h/3; } while (fabs((I0-I1)/I1)>eps); er=fabs((I0-I1)/I1); return I1; } //метод Гаусса double Gauss(double a, double b, double eps, int &n, double &er){ double t[6], A[6]; double sum=0; A[0]=0.17132450; A[1]=0.36076158; A[2]=0.46791394; A[5]=0.17132450; A[4]=0.36076158; A[3]=0.46791394; t[0]=-0.932469515; t[1]=-0.66120939; t[2]=-0.23861919; t[5]=0.932469515; t[4]=0.66120939; t[3]=0.23861919; double I0,I1=(f(a)+f(b))/2; for (int j=0;j<6;j++) { sum=0; for (int k=0;k<6;k++) { double x=(b+a*(2*j+1)+(b-a)*t[k])/(2*(j+1)); for (int i=0;i<=j;i++) sum+=A[k]*f(x+(b-a)*i/(j+1.0)); } I0=I1; I1=sum*(b-a)/(2*(j+1)); er=fabs((I0-I1)/I1); if (fabs((I0-I1)/I1) return I1; } int l=1; int r=1; if (fabs((I0-I1)/I1)>eps) I1=Gauss(a,a+(b-a)/2,eps,l,er)+Gauss(a+(b-a)/2,b,eps,r,er); er=fabs((I0-I1)/I1); n=l+r; return I1; } #pragma argsused int main(int argc, char* argv[]) { printf("Eps=?\n\r"); double eps; scanf("%lf",&eps); int n=0; double er; printf("Simpsons' method I=%f ",Simpson(0,2,eps,n,er)); printf("N=%i (er=%e)\n\r",n,er); n=1; printf("Gauss' method I=%f ",Gauss(0,2,eps,n,er)); printf("N=%i (er=%e)\n\r",n,er); getch(); return 0; } 1 2 |