数值分析 部分代码实现

新博客：https://wywwzjj.top/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define true  1
#define false 0
#define zero123 1.0E-20

double intf123(double x) { return  cos(x); } //用于积分
double fode123(double x, double y) { return -2 * x* y; }// 用于ode
double func123(double x){	return x*x*x + x * x - 3 * x - 3;}
double func1123(double x){return 3*x*x + 2 * x - 3;}

double erfen();
int Newton123(double x, double precision, int maxcyc);
double Trapez123(double a, double b, double eps);// 变步长复合梯形积分
double Romberg123(double a, double b, double eps);// Romberg积分
double Simpson123(double a, double b, double eps);// 变步长复合Simpson积分
void GSpivot123(double **A,double X[],int N); // 列选主元Gauss方法
void GSElim123(double **a,double b[],int n); //
double polyfit123(double A[], int m,double x0); // 拟合多项式函数
int GS123(double **a, double b[], int n, double x[], double eps); //Gauss-Seidel迭代法
int Jacobi123(double **a, double b[], int n, double x[], double eps);// Jacobi 迭代法
double Lagrange123(double x0, double x[], double y[], int n); //Lagrange 插值
double Newtontdd123(double x0, double x[], double Y[], int n);//Newton 插值
double mEuler123(double x, double y, double dt); //改进Euler方法
double RK4123(double x, double y, double dt); // 4阶Runge-Kutta方法

int main() 
{
	printf("选择程序类型-----\n");
	printf("0. 方程求根\n");
	printf("1. 方程组求解程序\n");
	printf("2. 积分程序\n");
	printf("3. 迭代法求解方程组程序\n");
	printf("4. 最小二乘法程序\n");
	printf("5. 插值程序\n");
	printf("6. 微分方程程序\n");
	printf("选择程序类型-----\n");
	int Select;
	int loopind = 1;

	while (loopind > 0) 
	{
		scanf("%d", &Select);
		if (Select == 0)//方程求根
		{
			double x, precision;
			int maxcyc;
			x = erfen();
			printf("输入最大迭代次数：");
			scanf("%d", &maxcyc);
			printf("迭代要求的精度：");
			scanf("%lf", &precision);
			if (Newton123(x, precision, maxcyc) == 1)//若函数返回值为1
			{
				printf("该值附近的根为：%lf\n", x);
			}
			else//若函数返回值为0
			{
				printf("迭代失败！\n");
			}
		}

		if (Select == 1)
		{
			int N = 3;
			double X[3]; //解向量
			double A1[3][4] = { {2, 2, -1, -4},	{1, -1, 0,0},{4, -2,-1, -6} };
			double **A;
			A = (double **)malloc(sizeof(double*) * N);
			for (int i = 0; i < N; i++)
				A[i] = (double *)malloc(sizeof(double) * (N+1));
			for (int i = 0; i < N; i++)
				for (int j = 0; j < (N+1); j++)
					A[i][j] = A1[i][j];
			printf("选择方法：1，列选主元； 2，Gauss消去法\n");
			int nSol;
			scanf("%d",&nSol);
			if (nSol==1) GSpivot123(A, X, N);
			else if (nSol == 2) GSElim123(A, X, N);
		}

		if (Select == 2) //积分程序
		{
			FILE *fpWrite = fopen("D:\\res.txt", "w");
			const double std = -1.120886483006822;
			double a0 = 1;
			double b0 = 6;
			double t1, t2, t3;
			double eps = 0.000001;
			t1 = Trapez123(a0, b0, eps);
			t2 = Simpson123(a0, b0, eps);
			t3 = Romberg123(a0, b0, eps);  // 不用分组
			fprintf(fpWrite, "变步长梯形积分结果：\n");
			for(int i = 0; i < 13; i++) {
				double n = pow(2.0, i);
				double per = (b0 - 1) / n;
				double sum1 = 0, sum2 = 0;
				for(int j = 0; j < (int)n; j++) {
					sum1 += Simpson123(1 + per * j, 1 + per * (j + 1), eps);
					sum2 += Trapez123(1 + per * j, 1 + per * (j + 1), eps);
				}
				fprintf(fpWrite, "k = %d\tSimpson积分结果: %6.6f\t误差：%e\tTrapez积分结果: %6.6f\t误差：%e\n",i, sum1, std-sum1, sum2, std-sum2);
			}
			fprintf(fpWrite, "\nRomberg积分结果：%6.6f\n", t3);
			fclose(fpWrite);
		}

		if (Select == 3) //迭代法求解方程组跟
		{
			int i, iter = 0;
			double eps;
			//double a[4][4] = { {7.0,2.0,1.0,-2.0},{9.0,15.0,3.0,-2.0},
			//		 {-2.0,-2.0,11.0,5.0},{1.0,3.0,2.0,13.0} };
			//double b[4] = { 4.0,7.0,-1.0,0.0 };
			FILE *fpA = fopen("D:/计算方法实习/jisuanshuzhi/main/A.txt", "r");
			if (fpA == NULL)
			{
				printf("未找到\n");
				return 0;
			}
			//读入数据
			double ndims = 9;
			double **IA;
			IA = (double **)malloc(sizeof(double*) * ndims);
			for (int i = 0; i < ndims; i++)
				IA[i] = (double *)malloc(sizeof(double) * ndims);
			for (int i = 0; i < ndims; i++)
				for (int j = 0; j < ndims; j++)
					fscanf(fpA, "%lf", &IA[i][j]);
			fclose(fpA);

			FILE *fpb = fopen("D:/计算方法实习/jisuanshuzhi/main/b.txt", "r");
			if (fpb == NULL) 
			{
			return 0;
			}
			double *bA;
			bA = (double *)malloc(sizeof(double) * ndims);
			for (int i = 0; i < ndims; i++)	fscanf(fpb, "%lf", &bA[i]);
			fclose(fpb);
			double *xiter; // xiteration
			xiter = (double *)malloc(sizeof(double) * ndims);
			eps = 0.000001;
			printf("迭代方法？ 1:GS123 ，2：Jacobi123\n");
			int itselect;
			scanf("%d", &itselect);

			if (itselect == 1)//高斯赛德尔
			{
				iter = GS123(IA, bA, ndims, xiter, eps);
				if (iter > 0)
				{
					printf("GS123法迭代次数：%d\n", iter);
					printf("GS123迭代法解为：\n");
					for (i = 0; i < ndims; i++) printf("x(%d)=%+6.6f\n", i, xiter[i]);
				}
			}
			else // Jacobi 方法
			{
				iter = Jacobi123(IA, bA, ndims, xiter, eps);
				if (iter > 0) 
				{
					printf("Jacobi123法迭代次数：%d\n", iter);
					printf("GS123迭代法解为：\n");
					for (i = 0; i < ndims; i++) printf("x(%d)=%+6.6f\n", i, xiter[i]);
				}
			}
		}

		if (Select == 4) //最小二乘法
		{
			double Xls[4] = { 1,2,3,4 }; // ls: least square
			double Yls[4] = { 1.95,2.4,2.83,3.3 };
			int N = 4;
			int m;
			printf("拟合次数(小于数据个数)：\n");
			scanf("%d", &m);
			double temp;
			double *S,*T,*Err,**C;
			double *LSA; // 拟合系数
			S = (double *)malloc(sizeof(double) * (2 * m+1));  // 动态分配数组
			T = (double *)malloc(sizeof(double) * (m + 1));

			Err = (double *)malloc(sizeof(double) * (m + 1));// 动态分配数组
			C = (double **)malloc(sizeof(double*) *(m+ 1));

			LSA = (double *)malloc(sizeof(double) * (m + 1));
			for (int i = 0; i <= m; i++)
				C[i] = (double *)malloc(sizeof(double) * (m+2));

			for (int i = 0; i <= 2 * m; i++)
			{
				temp = 0;
				for (int j = 0; j < N; j++)	temp += pow(Xls[j], i);
				S[i] = temp;
			}

			for (int i = 0; i <= m; i++)
			{
				temp = 0;
				for (int j = 0; j < N; j++)	temp += Yls[j] * pow(Xls[j], i);
				T[i] = temp;
			}
			for (int i = 0; i <= m; i++)
			{
				C[i][m + 1] = T[i];
				for (int j = 0; j <= m; j++) C[i][j] = S[i + j];
			}
			GSpivot123(C, LSA, m+1); // 利用列选主元消去法求解方程
			double err = 0;
			for (int i = 0; i < N; i++)
			{
				temp = polyfit123(LSA, m, Xls[i]);
				err += pow((Yls[i] - temp), 2);
			}
			err = sqrt(err / N);
			printf("%d次多项式拟合:\n", m);
			printf("\n最小均方误差%f:\n", err);
		}

		if (Select == 5)//插值
		{
			int n = 3;
			double ylag, ynt;
			double x0 = 50, X[3] = { 30,45,60 }, Y[3] = { 1 / 2.,sqrt(2.) / 2.,sqrt(3.) / 2. };
			ylag = Lagrange123(x0, X, Y, n);
			ynt = Newtontdd123(x0, X, Y, n);
			printf("在50°处的插值：\n");
			printf("\nLagrange 插值结果: %6.4f\n", ylag);
			printf("\nNewton 插值结果: %6.4f\n", ynt);
		}

		if (Select == 6)//微分方程
		{
			int nStep;
			double y0 = 1;
			double x0 = 0;
			// 初始值
			double y = y0;
			double x = x0;
			int odeselect ;
			printf("迭代方法？ 1:改进Euler ，2：RK4123\n");
			scanf("%d", &odeselect);
			printf("迭代步数:\n");
			scanf("%d", &nStep);
			double dt = 1. / nStep;
			printf("n=0: x=%6.3f y=%6.3f\n", x, y);
			if (odeselect == 1)
			{
				FILE *fpeuler = fopen("mEuer.out", "w");
				printf("修改Euler法:\n");
				for (int i = 1; i <= nStep; i++)
				{
					y = mEuler123(x, y, dt);
					x += dt;
					printf("n=%2d: x=%3.3f y=%8.6f\n", i, x, y);
					fprintf(fpeuler,"%2d %3.3f %12.8f\n", i, x, y);
				}
				fclose(fpeuler);
			}
			if (odeselect == 2)
			{
				printf("Runge-Kutta法\n");
				FILE *fprk = fopen("RK4123.out", "w");
				for (int i = 1; i <= nStep; i++)
				{
					y = RK4123(x, y, dt);
					x += dt;
					printf("n=%2d: x=%3.3f y=%8.6f\n", i, x, y);
					fprintf(fprk, "%2d %3.3f %12.8f\n", i, x, y);
				}
				fclose(fprk);
			}
		}
		printf("要继续吗？ 1:继续 ，0：退出");
		scanf("%d", &loopind);
		if (loopind == 1) printf("选择程序类型\n");
		else if (loopind == 0) system("pause");
	}
}
// 子函数
double erfen()
{
	double l, r, mid0123, ans0123, cnt;
	l = 1;
	r = 2;
	while (1) {//erfenfa
		mid0123 = (l + r) / 2;
		cnt = func123(mid0123);
		if (fabs(cnt) <= 0.0001) {
			ans0123 = mid0123;
			break;
		}
		if (cnt > 0)r = mid0123;
		else l = mid0123;
	}
	return ans0123;
}
int Newton123(double x, double precision, int maxcyc)//迭代次数  niudundiedai
{
	double x1, x0;
	int k;
	x0 = x;
	for (k = 0; k < maxcyc; k++)
	{
		if (func1123(x0) == 0.0)//若通过初值，函数返回值为0
		{
			printf("迭代过程中导数为0!\n");
			return 0;
		}
		x1 = x0 - func123(x0) / func1123(x0);//进行牛顿迭代计算
		if (fabs(x1 - x0) < precision || fabs(func123(x1)) < precision)//达到结束条件
		{
			x = x1;//返回结果
			return 1;
		}
		else//未达到结束条件
		{
			x0 = x1;//准备下一次迭代
		}
	}
	printf("迭代次数超过预期！\n");//迭代次数达到，仍没有达到精度
	return 0;
}
double Trapez123(double a, double b, double eps)//变步长复合Simpson积分
{
	int n, k;
	double fa, fb, h, t1, p, s, x, t;
	fa = intf123(a);
	fb = intf123(b);
	n = 1; h = b - a;
	t1 = h * (fa + fb) / 2.0;
	p = eps + 1.0;
	while (p >= eps)
	{
		s = 0.0;
		for (k = 0; k <= n - 1; k++)
		{
			x = a + (k + 0.5)*h;
			s = s + intf123(x);
		}
		t = (t1 + h * s) / 2.0;
		p = fabs(t1 - t);
		t1 = t; n = n + n; h = h / 2.0;
	}
	return(t);
}
double Romberg123(double a, double b, double eps)//longbeige
{
	int m, n, i, k;
	double y[10], h, ep, p, x, s, q;
	h = b - a;
	y[0] = h * (intf123(a) + intf123(b)) / 2.0;
	m = 1; n = 1; ep = eps + 1.0;
	while ((ep >= eps) && (m <= 9))
	{
		p = 0.0;
		for (i = 0; i <= n - 1; i++)
		{
			x = a + (i + 0.5)*h;
			p = p + intf123(x);
		}
		p = (y[0] + h * p) / 2.0;
		s = 1.0;
		for (k = 1; k <= m; k++)
		{
			s = 4.0*s;
			q = (s*p - y[k - 1]) / (s - 1.0);
			y[k - 1] = p; p = q;
		}
		ep = fabs(q - y[m - 1]);
		m = m + 1; y[m - 1] = q;
		n = n + n; h = h / 2.0;
	}
	return(q);
}
double Simpson123(double a, double b, double eps)//xinpusheng
{
	int n, k;
	double h, t1, t2, s1, s2, ep, p, x;
	n = 1;
	h = b - a;
	t1 = h * (intf123(a) + intf123(b)) / 2.0;
	s1 = t1;
	ep = eps + 1.0;
	while (ep >= eps)
	{
		p = 0.0;
		for (k = 0; k <= n - 1; k++)
		{
			x = a + (k + 0.5)*h;
			p = p + intf123(x);
		}
		t2 = (t1 + h * p) / 2.0;
		s2 = (4.0*t2 - t1) / 3.0;
		ep = fabs(s2 - s1);
		t1 = t2; s1 = s2;
		n = n + n;
		h = h / 2.0;
	}
	return(s2);
}

double polyfit123(double A[], int m,double x0) // 拟合多项式函数
{
	double temp = 0;
	for (int i = 0; i <= m; i++)
		temp += A[i]*pow(x0, i);
	return temp;
}
void GSpivot123(double **A,double X[],int N)
{
	int M, ICHG, I, NN, IMAX, J, JJ, IP, JP, NCOPY, I1, J1, N1, K, N2, KK, OK = true;
	int *NROW;
	double AMAX, XM, SUM;  /* AMAX save  a_{pk} = max_{k<=i<=n} |a_{ik}^{(k)}| */
	//int NROW[N];
	NROW = (int *)malloc(sizeof(int) * N);// 保存主元变换顺序

	M = N + 1;
	for (I = 1; I <= N; I++) NROW[I - 1] = I;
	NN = N - 1;
	ICHG = 0; // index change
	I = 1;
	/* STEP 2 */
	while ((OK) && (I <= NN)) {
		/* STEP 3 */
		// For the scaled partial pivoting, AMAX should be scaled by the scale factor
		IMAX = NROW[I - 1];
		AMAX = fabs(A[IMAX - 1][I - 1]);
		IMAX = I;
		JJ = I + 1;
		for (IP = JJ; IP <= N; IP++)
		{
			JP = NROW[IP - 1];
			if (fabs(A[JP - 1][I - 1]) > AMAX)
			{
				AMAX = fabs(A[JP - 1][I - 1]);
				IMAX = IP;
			}

		}
		if (AMAX <= zero123) OK = false;
		else {
			if (NROW[I - 1] != NROW[IMAX - 1])
			{
				ICHG = ICHG + 1;
				NCOPY = NROW[I - 1];
				NROW[I - 1] = NROW[IMAX - 1];
				NROW[IMAX - 1] = NCOPY;
			}
			I1 = NROW[I - 1];

			for (J = JJ; J <= N; J++)
			{
				J1 = NROW[J - 1];
				XM = A[J1 - 1][I - 1] / A[I1 - 1][I - 1];
				for (K = JJ; K <= M; K++)
					A[J1 - 1][K - 1] = A[J1 - 1][K - 1] - XM * A[I1 - 1][K - 1];
				A[J1 - 1][I - 1] = 0.0;
			}
		}
		I++;
	}
	if (OK)
	{
		N1 = NROW[N - 1];
		if (fabs(A[N1 - 1][N - 1]) <= zero123) OK = false;
		else
		{
			// 回代过程
			X[N - 1] = A[N1 - 1][M - 1] / A[N1 - 1][N - 1];
			for (K = 1; K <= NN; K++) {
				I = NN - K + 1;
				JJ = I + 1;
				N2 = NROW[I - 1];
				SUM = 0.0;
				for (KK = JJ; KK <= N; KK++) {
					SUM = SUM - A[N2 - 1][KK - 1] * X[KK - 1];
				}
				X[I - 1] = (A[N2 - 1][N] + SUM) / A[N2 - 1][I - 1];
			}
			// 输出解
			printf("解向量：\n");
			for (I = 1; I <= N; I++)
			{
				printf("  %12.8f", X[I - 1]);
			}
			printf("\n");
		}
	}
	if (!OK) printf("系统没有解\n");
}
void GSElim123(double **a,double b[],int n)
{
	int i,j,k;
	double c,d;
	printf("\n Gauss消去法求解\n");

	// apply elimination
	for(k=0;k<n-1;k++)
	{
		for(i=k+1;i<n;i++)
		{
			d=a[i][k];
			for(j=0;j<n+1;j++)
			{
				c=(a[k][j]/a[k][k]);
				a[i][j]-=c*d;
			}
		}
	}
	printf("\n 消除后的矩阵:\n");
	for(i=0;i<n;i++)
	{
		printf("\n");
		for(j=0;j<n+1;j++)
			printf("%8.6f\t",a[i][j]);
	}
	printf("\n");

	for(i=n-1;i>=0;i--)
	{
		b[i]=a[i][n];
		d=0;
		for(k=1;k<n-i;k++)
			d+=(a[i][n-k]*b[n-k]);
		b[i]=(b[i]-d)/a[i][i];
	}
	printf("解向量:\n  ");
	for(i=0;i<n;i++)
		printf("  %12.8f", b[i]);
	printf("\n");
}
int GS123(double **a, double b[], int n, double x[], double eps)
{
	int i, j,iter=0;

	double p, t, s, q;
	for (i = 0; i <= n - 1; i++)
	{
		p = 0.0;
		x[i] = 0.0;
		for (j = 0; j <= n - 1; j++)
			if (i != j)
			{
				p = p + fabs(a[i][j]);
			}
		if (p >= fabs(a[i][i]))
		{
			printf("fail\n");
			return(-1);
		}
	}
	p = eps + 1.0;
	while (p >= eps)
	{
		p = 0.0;
		for (i = 0; i <= n - 1; i++)
		{
			t = x[i];
			s = 0.0;
			for (j = 0; j <= n - 1; j++)
			{
				if (j != i) s = s + a[i][j] * x[j];
			}

			x[i] = (b[i] - s) / a[i][i];
			q = fabs(x[i] - t) / (1.0 + fabs(x[i]));
			if (q > p) p = q;
		}
		iter = iter + 1;
	}
	return iter;
}
int Jacobi123(double **a, double b[], int n, double x[], double eps)
{
	int i, j, iter = 0;
	double p, t, s, q;
	double *x0;
	x0 = (double *)malloc(sizeof(double) * n);
	//double **a;

	for (i = 0; i <= n - 1; i++)
	{
		p = 0.0;
		x[i] = 0.0;
		for (j = 0; j <= n - 1; j++)
			if (i != j)
			{
				p = p + fabs(a[i][j]);
			}
		if (p >= fabs(a[i][i]))
		{
			printf("fail\n");
			return(-1);
		}
	}
	p = eps + 1.0;
	while (p >= eps)
	{
		p = 0.0;
		for (i = 0; i <= n - 1; i++)
		{
			x0[i] = x[i];
		}
		for (i = 0; i <= n - 1; i++)
		{
			t = x[i];
			s = 0.0;
			for (j = 0; j <= n - 1; j++)
			{
				if (j != i) s = s + a[i][j] * x0[j];
			}

			x[i] = (b[i] - s) / a[i][i];
			q = fabs(x[i] - t) / (1.0 + fabs(x[i]));
			if (q > p) p = q;
		}
		iter = iter + 1;
	}
	return iter;
}

double Lagrange123(double x0, double x[], double y[], int n)
{
	double sum = 0;
	for (int i = 0; i < n; i++) {
		double temp = 1.0;
		for (int j = 0; j < n; j++)
		{
			if (j == i)	continue;
			temp *= (x0 - x[j]) / (x[i] - x[j]);
		}
		sum += y[i] * temp;
	}
	return sum;
}
double Newtontdd123(double x0, double x[], double Y[], int n)
{
	int i;
	double **y;
	y = (double **)malloc(sizeof(double*) * n);
	for (int i = 0; i < n; i++)
		y[i] = (double *)malloc(sizeof(double) * n);
	//double y[][n];
	for (i = 0; i < n; i++)
		y[0][i] = Y[i];
	double sum = y[0][0];                         //插值计算
	double x_xi = x0 - x[0];
	for (int i = 1; i < n; i++) {
		for (int j = 0; j < n - i; j++) {
			double a1 = x[j];                     //使用递推方法求解插值
			double a2 = x[i + j];                 //i+j获取副对角线元素下标
			y[i][j] = (y[i - 1][j + 1] - y[i - 1][j]) / (a2 - a1);
		}
		sum += y[i][0] * x_xi;
		x_xi *= (x0 - x[i]);
	}
	return sum;
}

double mEuler123(double x, double y, double dt)
{
	double predictor = y + (dt * fode123(x, y));
	double k1 = fode123(x, y);
	double k2 = fode123(x + dt, predictor);
	return y += dt * 0.5 * (k1 + k2);
}
double RK4123(double x, double y, double dt)
{
	double k1 = fode123(x, y);
	double k2 = fode123(x + dt / 2, y + dt / 2 * k1);
	double k3 = fode123(x + dt / 2, y + dt / 2 * k2);
	double k4 = fode123(x + dt, y + dt * k3);
	return y += dt / 6 * (k1 + 2 * k2 + 2 * k3 + k4);
}