工程计算实验代码:

最新推荐文章于 2022-08-28 18:34:32 发布
最新推荐文章于 2022-08-28 18:34:32 发布
阅读量488 收藏 2
点赞数 1
分类专栏： 工程计算
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qq_42683011/article/details/105724065
版权

实验一:

实验二:

#include <iostream>

using namespace std;

//调试用, 增加更多的输出 
//#define DEBUG

void Output1D(double matrix[], int row);

int main() {
	int n = 0;
	while(n<15) {
		n++;
		printf(
		    "-------------------------------------------------\n"
		    "第%d次循环:  n=%d\n"
		    "\n",n,n);
//		scanf("%d", &n);
		double H[n][n] = { 0 };
		double b[n]= {0};
		double l[n][n]= {0};
		double d[n]= {0};
		double y[n]= {0};
		double x[n]= {0};

		//初始化希尔伯特矩阵H
		for (int i = 0; i < n; ++i) {
			for (int j = 0; j < n; ++j) {
				H[i][j] = 1.0 / (i + 1 + j + 1 - 1);
			}
		}
#ifdef DEBUG
		//输出h:
		printf("H:\n");
		for(int i=0; i<n; ++i) {
			Output1D(H[i],n);
			putchar('\n');
		}
		putchar('\n');
		putchar('\n');
#endif
		//求得矩阵b
		for (int i = 0; i < n; ++i) {
			for (int j = 0; j < n; ++j) {
				b[i] += H[i][j];
			}
		}
#ifdef DEBUG
		//输出b
		printf("b:\n");
		Output1D(b,n);
		putchar('\n');
		putchar('\n');
#endif

		//初始化矩阵l
		for(int i=0; i<n; i++) {
			for(int j=0; j<n; j++) {
				l[i][j]=0;
			}
		}
		for (int i = 0; i < n; ++i) {
			l[i][i]=1;
		}
#ifdef DEBUG
		//输出l:
		printf("L:\n");
		for(int i=0; i<n; ++i) {
			Output1D(l[i],n);
			putchar('\n');
		}
		putchar('\n');
#endif

		//对希尔伯特矩阵H进行分解
		//H=LDL(T)
		int row =n;
		for (int k = 0; k < row; ++k) {
			d[k] = H[k][k];
			for (int i = 0; i < k; ++i) {
				d[k] -= l[k][i] * l[k][i] * d[i];
			}
			for (int i = k + 1; i < row; ++i) {
				double temp = H[k][i];
				for (int j = 0; j < k; ++j) {
					temp -= l[k][j] * d[j] * l[i][j];
				}
				l[i][k] = temp / (l[k][k] * d[k]);
			}
		}
#ifdef DEBUG
		//输出d:
		printf("D:\n");
		Output1D(d,n);
		putchar('\n');
		putchar('\n');
#endif

#ifdef DEBUG
		//输出l
		printf("L:\n");
		for(int i=0; i<n; ++i) {
			Output1D(l[i],n);
			putchar('\n');
		}
		putchar('\n');
#endif

		//对向量y进行求解
		//Ly=b
		y[0] = b[0];
		for (int i = 1; i < row; ++i) {
			y[i] = b[i];
			for (int k = 0; k < i; ++k) {
				y[i] -= l[i][k] * y[k];
			}
		}
		//输出y:
		printf("Y:\n");
		Output1D(y,n);
		putchar('\n');
		putchar('\n');

		//对向量x进行求解
		//L(t)x=y
		x[row - 1] = y[row - 1] / d[row - 1];
		for (int i = row - 2; i >= 0; --i) {
			x[i] = y[i] / d[i];
			for (int k = i + 1; k < row; ++k) {
				x[i] -= l[k][i] * x[k];
			}
		}
		//输出x:
		printf("X:\n");
		Output1D(x,n);
		putchar('\n');
		putchar('\n');

		//计算剩余向量
		//r=b-Hx
		double r[n]= {0};
		double temp=0;
		for(int i=0; i<n; i++) {
			temp=0;
			for(int j=0; j<n; ++j) {
				temp+=H[i][j]*x[j];
			}
			r[i]=b[i]-temp;
		}
		//输出r:
		printf("r:\n");
		Output1D(r,n);
		putchar('\n');
		putchar('\n');

		//计算与x=(1,1, ... ,1)T 的差值
		//diff_x=x'-x
		double diff_x[n]= {0};
		for(int i=0; i<n; ++i) {
			diff_x[i]=x[i]-1;
		}
		//输出diff_x:
		printf("diff_x:\n");
		Output1D(diff_x,n);
		putchar('\n');
		putchar('\n');

	}
	return 0;
}

void Output1D(double matrix[], int n) {
	for (int i = 0; i < n; ++i) {
		printf("%.8lf ", matrix[i]);
	}
	return;
}

实验三:

源码:

Jacobi:

#include <iostream>
#include <cmath>
#include <cstring>

//#define DEBUG		//DEBUG用, 增加更多的输出
#define PRO_1		//第一题时需要定义


using namespace std;

const int N=30;
const int INF=10000;

void iniA(double a[N][N], int n) ;
void iniB(double b[N],int n);
void input1D(double vec[N], int n) ;
void input2D(double matrix[N][N],int n);
void output1D(double vec[N],int n);
void output2D(double matrix[N][N], int n);
void VectorMinus(double a[N],double b[N],double ans[N],int n);
void MatrixMtimes(double a[N][N], double b[N][N],double ans[N][N],int n) ;
void VectorCopy(double src[N],double dst[N],int n);
void MatrixCopy(double src[N][N],double dst[N][N],int n) ;
double Norm2(double a[N],int n) ;


int main(int argc, char* argv[]) {

	double ma0[N][N]= {0};
	double coefficientB[N]= {0};
	double x_k1[N]= {0};
	double x_k0[N]= {0};
	double x_diff[N]= {0};
	double eps=0.0;

//------------------------------------
//第一题:
#ifdef PRO_1
	//默认矩阵为方阵
	int row=0;
	printf(
	    "输入矩阵的行列:\n"
	    "row = column =  ");
	scanf("%d",&row);

	printf("输入矩阵A:\n");
	input2D(ma0,row);

	printf("输入系数矩阵B:\n");
	input1D(coefficientB,row);


	printf("输入初始X0 = ");
	for(int i=0; i<row; ++i) {
		scanf("%lf",&x_k0[i]);
	}

	printf("输入终止条件: eps = ");
	scanf("%lf",&eps);

	int k=0;
	bool isDiverge=0;
	do {

#ifdef DEBUG
		printf(
		    "第%d次循环\n"
		    "X_k0=\n"
		    ,k);
		output1D(x_k0,row);
#endif
		for(int i=0; i<row; ++i) {
			double sum=0.0;
			for(int j=0; j<row; ++j) {
				sum+=ma0[i][j]*x_k0[j];
			}
			sum-=ma0[i][i]*x_k0[i];
			sum-=coefficientB[i];
			sum/=(-ma0[i][i]);
			x_k1[i]=sum;
		}
		double tempAns[N]= {0};
		VectorMinus(x_k1,x_k0,tempAns,row);
		VectorCopy(tempAns,x_diff,row);
		VectorCopy(x_k1,x_k0,row);
#ifdef DEBUG
		printf("x_diff=\n");
		output1D(x_diff,row);
		printf(
		    "x_diff的2范数 = \n"
		    "%14.10lf\n"
		    ,Norm2(x_diff,row));
#endif

		//当Jacobi的结果超过既定的最大值INF时, 判定为Jacobi发散:
		for(int i=0; i<row; ++i) {
			if(fabs(x_k1[i])>=INF) {
				isDiverge=true;
				break;
			}
		}
		if(isDiverge) {
			break;
		}
	} while(Norm2(x_diff,row)>=eps);


	if(isDiverge) {
		isDiverge=false;
		printf("Jacobi发散....\n");
	} else {
		printf("最终结果: x_k1=\n");
		output1D(x_k1,row);
	}

///--------------------------------------------------
//第二题:
#else
	int ROW=0, row=0;
	printf("输入方阵的最大行数: ROW = ");
	scanf("%d",&ROW);
	while(row<ROW) {
		++row;
		printf("row = %d\n",row);
		iniA(ma0,row);
		iniB(coefficientB,row);
		eps=0.000001;

#ifdef DEBUG
		printf("矩阵A:\n");
		output2D(ma0,row);
		printf("系数矩阵B:\n");
		output1D(coefficientB,row);
#endif

		int k=0;
		const int DIVERGE=100;
		bool isDiverge=0;
		do {
			//当迭代次数>DIVERGE时判定为发散
			if(++k>DIVERGE) {
				isDiverge=true;
				break;
			}

#ifdef DEBUG
			printf(
			    "第%d次循环\n"
			    "X_k0=\n"
			    ,k);
			output1D(x_k0,row);
#endif
			for(int i=0; i<row; ++i) {
				double sum=0.0;
				for(int j=0; j<row; ++j) {
					sum+=ma0[i][j]*x_k0[j];
				}
				sum-=ma0[i][i]*x_k0[i];
				sum-=coefficientB[i];
				sum/=(-ma0[i][i]);
				x_k1[i]=sum;
			}
			double tempAns[N]= {0};
			VectorMinus(x_k1,x_k0,tempAns,row);
			VectorCopy(tempAns,x_diff,row);
			VectorCopy(x_k1,x_k0,row);
#ifdef DEBUG
			printf("x_diff=\n");
			output1D(x_diff,row);
			printf(
			    "x_diff的2范数 = \n"
			    "%14.10lf\n"
			    ,Norm2(x_diff,row));
#endif


			//当Jacobi的结果超过既定的最大值INF时, 判定为Jacobi发散:
			for(int i=0; i<row; ++i) {
				if(fabs(x_k1[i])>=INF) {
					isDiverge=true;
					break;
				}
			}
			if(isDiverge) {
				break;
			}
		} while(Norm2(x_diff,row)>=eps);

		if(isDiverge) {
			isDiverge=false;
			printf("Jacobi发散....\n");
		} else {
			printf(
			    "最终结果: \n"
			    "经过%d次迭代\n"
			    "最终结果: x_k1=\n"
			    ,k);
			output1D(x_k1,row);
		}
	}

#endif

	return 0;
}

//向量减法:
void VectorMinus(double a[N],double b[N],double ans[N],int n) {
	const double MINIMUM=pow(10,-11);
	for(int i=0; i<n; ++i) {
		ans[i]=a[i]-b[i];
		if(ans[i]<MINIMUM) {
			ans[i]=0;
		}
	}
	return ;
}

//方阵乘法:
void MatrixMtimes(double a[N][N], double b[N][N],double ans[N][N],int n) {
	const double MINIMUM=pow(10,-11);
	memset(ans,0,sizeof(double)*N*N);
	for(int i=0; i<n; ++i) {
		for(int j=0; j<n; ++j) {
			for(int k=0; k<n; ++k) {
				ans[i][j]+=a[i][k]*b[k][j];
			}
			//当值小于10^-11时, 视为=0
			if(ans[i][j]<=MINIMUM) {
				ans[i][j]=0;
			}
		}
	}
	return ;
}

//向量复制:
void VectorCopy(double src[N],double dst[N],int n) {
	for(int i=0; i<n; ++i) {
		dst[i]=src[i];
	}
	return ;
}


//求向量的2范数:
double Norm2(double a[N],int n) {
	double sum=0.0;
	for(int i=0; i<n; ++i) {
		sum+=fabs(a[i]*a[i]);
	}
	return sqrt(sum);
}


void iniA(double a[N][N], int n) {
	for(int i=0; i<n; ++i) {
		for(int j=0; j<n; ++j) {
			a[i][j]= 1.0/(i+j+1);
		}
	}
	return ;
}

void iniB(double b[N],int n) {
	double sum=0;
	for(int i=0; i<n; ++i) {
		sum=0;
		for(int j=0; j<n; ++j) {
			sum+=1.0/(i+j+1);
		}
		b[i]=sum;
	}
	return ;
}

void input1D(double vec[N], int n) {
	for(int i=0; i<n; ++i) {
		scanf("%lf",&vec[i]);
	}
	return ;
}

void input2D(double matrix[N][N],int n) {
	for(int i=0; i<n; ++i) {
		input1D(matrix[i],n);
	}
	return ;
}

void output1D(double vec[N],int n) {
	for(int i=0; i<n; ++i) {
		printf("%+18.10lf ",vec[i]);
	}
	putchar('\n');
	return ;
}

void output2D(double matrix[N][N], int n) {
	for(int i=0; i<n; ++i) {
		output1D(matrix[i],n);
	}
	putchar('\n');
	return ;
}

SOR:

#include <iostream>
#include <cmath>
#include <cstring>

//#define DEBUG		//DEBUG用, 增加更多的输出
//#define PRO_1		//第一题时需要定义


using namespace std;

const int N=30;
const int INF=10000;

void iniA(double a[N][N], int n) ;
void iniB(double b[N],int n);
void input1D(double vec[N], int n) ;
void input2D(double matrix[N][N],int n);
void output1D(double vec[N],int n);
void output2D(double matrix[N][N], int n);
void VectorMinus(double a[N],double b[N],double ans[N],int n);
void VectorCopy(double src[N],double dst[N],int n);
double Norm2(double a[N],int n) ;

int main(void) {

	double ma0[N][N]= {0};
	double coefficientB[N]= {0};
	double x_k1[N]= {0};
	double x_k0[N]= {0};
	double x_diff[N]= {0};
	double eps=0.0;
	double omega=0.0;

//----------------------------------------
//第一题:
#ifdef PRO_1
	//默认矩阵为方阵
	int row=0;
	printf(
	    "输入矩阵的行列:\n"
	    "row = column =  ");
	scanf("%d",&row);

	printf("输入矩阵A:\n");
	input2D(ma0,row);

	printf("输入系数矩阵B:\n");
	input1D(coefficientB,row);


	printf("输入初始X0 = ");
	for(int i=0; i<row; ++i) {
		scanf("%lf",&x_k0[i]);
	}

	printf("输入终止条件: eps = ");
	scanf("%lf",&eps);

	printf("输入松弛因子: omega = ");
	scanf("%lf",&omega);

	int k=0;
	bool isDiverge=0;

	do {

#ifdef DEBUG
		printf(
		    "第%d次循环\n"
		    "X_k0=\n"
		    ,k);
		output1D(x_k0,row);
#endif

		for(int i=0; i<row; ++i) {
			double sum=0.0;
			for(int j=0; j<i; ++j) {
				sum-=ma0[i][j]*x_k1[j];
			}
			for(int j=i; j<row; ++j) {
				sum-=ma0[i][j]*x_k0[j];
			}
			sum+=coefficientB[i];
			sum*=omega;
			sum/=ma0[i][i];
			sum+=x_k0[i];
			x_k1[i]=sum;
		}
		VectorMinus(x_k1,x_k0,x_diff,row);
		VectorCopy(x_k1,x_k0,row);
#ifdef DEBUG
		printf("x_k1=\n");
		output1D(x_k1,row);
		printf("x_diff=\n");
		output1D(x_diff,row);
		printf(
		    "x_diff的2范数 = \n"
		    "%14.10lf\n"
		    ,Norm2(x_diff,row));
#endif

		//当SOR的结果超过既定的最大值INF时, 判定为发散:
		for(int i=0; i<row; ++i) {
			if(fabs(x_k1[i])>=INF) {
				isDiverge=true;
				break;
			}
		}
		if(isDiverge) {
			break;
		}
	} while(Norm2(x_diff,row)>=eps);
	if(isDiverge) {
		isDiverge=false;
		printf("SOR发散....\n");
	} else {
		printf(
		    "最终结果: \n"
		    "经过%d次迭代\n"
		    "最终结果: x_k1=\n"
		    ,k);
		output1D(x_k1,row);
	}
//----------------------------------------
//第二题:
#else
	int ROW=0, row=0;
	printf("输入方阵的最大行数: ROW = ");
	scanf("%d",&ROW);
	printf("输入松弛因子: omega = ");
	scanf("%lf",&omega);
	while(row<ROW) {
		++row;
		printf("row = %d\n",row);
		iniA(ma0,row);
		iniB(coefficientB,row);
		eps=0.000001;

#ifdef DEBUG
		printf("矩阵A:\n");
		output2D(ma0,row);
		printf("系数矩阵B:\n");
		output1D(coefficientB,row);
#endif

		int k=0;
		bool isDiverge=false;

		do {
			k++;


#ifdef DEBUG
			printf(
			    "第%d次循环\n"
			    "X_k0=\n"
			    ,k);
			output1D(x_k0,row);
#endif

			for(int i=0; i<row; ++i) {
				double sum=0.0;
				for(int j=0; j<i; ++j) {
					sum-=ma0[i][j]*x_k1[j];
				}
				for(int j=i; j<row; ++j) {
					sum-=ma0[i][j]*x_k0[j];
				}
				sum+=coefficientB[i];
				sum*=omega;
				sum/=ma0[i][i];
				sum+=x_k0[i];
				x_k1[i]=sum;
			}
			VectorMinus(x_k1,x_k0,x_diff,row);
			VectorCopy(x_k1,x_k0,row);
#ifdef DEBUG
			printf("x_k1=\n");
			output1D(x_k1,row);
			printf("x_diff=\n");
			output1D(x_diff,row);
			printf(
			    "x_diff的2范数 = \n"
			    "%14.10lf\n"
			    ,Norm2(x_diff,row));
#endif

			//当SOR的结果超过既定的最大值INF时, 判定为发散:
			for(int i=0; i<row; ++i) {
			if(fabs(x_k1[i])>=INF) {
					isDiverge=true;
					break;
				}
			}
			if(isDiverge) {
				break;
			}
		} while(Norm2(x_diff,row)>=eps);

		if(isDiverge) {
			isDiverge=false;
			printf("SOR发散....\n");
		} else {
			printf(
			    "最终结果: \n"
			    "经过%d次迭代\n"
			    "最终结果: x_k1=\n"
			    ,k);
			output1D(x_k1,row);
			printf(
			    "\n"
			    "-------------------------------------------\n");
		}

	}
#endif
	return 0;
}

//向量减法:
void VectorMinus(double a[N],double b[N],double ans[N],int n) {
	const double MINIMUM=pow(10,-11);
	for(int i=0; i<n; ++i) {
		ans[i]=a[i]-b[i];
		if(ans[i]<MINIMUM) {
			ans[i]=0;
		}
	}
	return ;
}

//向量复制:
void VectorCopy(double src[N],double dst[N],int n) {
	for(int i=0; i<n; ++i) {
		dst[i]=src[i];
	}
	return ;
}

//求向量的2范数:
double Norm2(double a[N],int n) {
	double sum=0.0;
	for(int i=0; i<n; ++i) {
		sum+=fabs(a[i]*a[i]);
	}
	return sqrt(sum);
}

void input1D(double vec[N], int n) {
	for(int i=0; i<n; ++i) {
		scanf("%lf",&vec[i]);
	}
	return ;
}

void input2D(double matrix[N][N],int n) {
	for(int i=0; i<n; ++i) {
		input1D(matrix[i],n);
	}
	return ;
}

void output1D(double vec[N],int n) {
	for(int i=0; i<n; ++i) {
		printf("%+18.10lf ",vec[i]);
	}
	putchar('\n');
	return ;
}

void output2D(double matrix[N][N], int n) {
	for(int i=0; i<n; ++i) {
		output1D(matrix[i],n);
	}
	putchar('\n');
	return ;
}

void iniA(double a[N][N], int n) {
	for(int i=0; i<n; ++i) {
		for(int j=0; j<n; ++j) {
			a[i][j]= 1.0/(i+j+1);
		}
	}
	return ;
}

void iniB(double b[N],int n) {
	double sum=0;
	for(int i=0; i<n; ++i) {
		sum=0;
		for(int j=0; j<n; ++j) {
			sum+=1.0/(i+j+1);
		}
		b[i]=sum;
	}
	return ;
}

CG:

#include <iostream>
#include <cmath>
#include <cstring>

//#define DEBUG		//DEBUG用, 增加更多的输出
//#define PRO_1		//第一题时需要定义


using namespace std;

const int N=30;
const int INF=10000;

void iniA(double a[N][N], int n) ;
void iniB(double b[N],int n);
void input1D(double vec[N], int n) ;
void input2D(double matrix[N][N],int n);
void output1D(double vec[N],int n);
void output2D(double matrix[N][N], int n);
double VectorMtimes(double a[N], double b[N],int n);
void VectorPlus(double a[N],double b[N], double ans[N], int n);
void VectorMinus(double a[N],double b[N],double ans[N],int n) ;
void MatrixMtimes(double a[N][N], double b[N],double ans[N],int n) ;
void VectorCopy(double src[N],double dst[N],int n) ;
double Norm2(double a[N],int n) ;
double NormInf(double a[N],int n);



int main(void) {

	double a[N][N]= {0};
	double b[N]= {0};
	double x_k1[N]= {0};
	double x_k0[N]= {0};
	double eps=0.0;
	double r_k0[N]= {0};
	double p_k0[N]= {0};
	double beta=0.0;
	double alpha=0.0;
	double tempAns[N]= {0};
	double tempAns2[N]= {0};


//--------------------------------------
//第一题:
#ifdef PRO_1
	//默认矩阵为方阵
	int row=0;
	printf(
	    "输入矩阵的行列:\n"
	    "row = column =  ");
	scanf("%d",&row);

	printf("输入矩阵A:\n");
	input2D(a,row);

	printf("输入系数矩阵B:\n");
	input1D(b,row);

#ifdef DEBUG
	printf("输入的矩阵A为:\n");
	output2D(a,row);
	printf("输入的系数b为:\n");
	output1D(b,row);
#endif

	printf("输入初始X0 = ");
	for(int i=0; i<row; ++i) {
		scanf("%lf",&x_k0[i]);
	}

	printf("输入终止条件: eps = ");
	scanf("%lf",&eps);

//计算初次迭代的各个初值:

	//r=b-Ax
	memset(tempAns,0,sizeof(double)*N);		//中间结果tempAns清零
	MatrixMtimes(a,x_k0,tempAns,row);
	VectorMinus(b,tempAns,r_k0,row);

	//p=r
	VectorCopy(r_k0,p_k0,row);

	//α= r^T*r/(p^T*A*p)
	memset(tempAns,0,sizeof(double)*N);		//中间结果tempAns清零
	MatrixMtimes(a,p_k0,tempAns,row);
	alpha=VectorMtimes(r_k0,r_k0,row)/VectorMtimes(p_k0,tempAns,row);

	//x'=x+αp
	memset(tempAns,0,sizeof(double)*N);		//中间结果tempAns清零
	for(int i=0; i<row; ++i) {
		tempAns[i]=p_k0[i]*alpha;
	}
	VectorPlus(x_k0,tempAns,x_k1,row);

#ifdef DEBUG
	printf(
	    "初始值:\n"
	    "r=\n"
	);
	output1D(r_k0,row);
	printf("p = \n");
	output1D(p_k0,row);
	printf("x_k0 = \n");
	output1D(x_k0,row);
	printf("x_k1 = \n");
	output1D(x_k1,row);
	printf("alpha = %lf\n",alpha);
#endif

	int k=0;
	double diff=0.0;
	double isDiverge=0;

	do {
		k++;

		//r=b-Ax
		memset(tempAns,0,sizeof(double)*N);		//中间结果tempAns清零
		MatrixMtimes(a,x_k1,tempAns,row);
		VectorMinus(b,tempAns,r_k0,row);
#ifdef DEBUG
		printf("tempAns = \n");
		output1D(tempAns,row);
#endif

		//β=r^TAp/p^TAp
		memset(tempAns,0,sizeof(double)*N);		//中间结果tempAns清零
		MatrixMtimes(a,p_k0,tempAns,row);
		beta=-VectorMtimes(r_k0,tempAns,row)/VectorMtimes(p_k0,tempAns,row);

		//p'=r+βp
		memset(tempAns,0,sizeof(double)*N);		//中间结果tempAns清零
		for(int i=0; i<row; ++i) {
			tempAns[i]=p_k0[i]*beta;
		}
		VectorPlus(r_k0,tempAns,p_k0,row);

		//α=r^Tp/p^TAp
		memset(tempAns,0,sizeof(double)*N);		//中间结果tempAns清零
		MatrixMtimes(a,p_k0,tempAns,row);
		alpha=VectorMtimes(r_k0,p_k0,row)/VectorMtimes(p_k0,tempAns,row);

		//x'=x+αp
		memset(tempAns,0,sizeof(double)*N);		//中间结果tempAns清零
		memset(tempAns2,0,sizeof(double)*N);	//中间结果tempAns2清零
		for(int i=0; i<row; ++i) {
			tempAns[i]=p_k0[i]*alpha;
		}
		VectorPlus(x_k1,tempAns,tempAns2,row);
		VectorCopy(tempAns2,x_k1,row);

		//使用diff残差判定:
		memset(tempAns,0,sizeof(double)*N);		//中间结果tempAns清零
		memset(tempAns2,0,sizeof(double)*N);	//中间结果tempAns2清零
		MatrixMtimes(a,x_k1,tempAns,row);
		VectorMinus(b,tempAns,tempAns2,row);
		diff=NormInf(tempAns2,row)/NormInf(b,row);

#ifdef DEBUG
		printf("第%d次循环\n",k);
		printf("r = \n");
		output1D(r_k0,row);
		printf("p = \n");
		output1D(p_k0,row);
		printf("x_k0 = \n");
		output1D(x_k0,row);
		printf("x_k1 = \n");
		output1D(x_k1,row);
		printf("alpha = %lf\n",alpha);
		printf("beta = %lf\n",beta);
		printf("diff = %lf\n",diff);
#endif
		//当CG的结果超过既定的最大值INF时, 判定为发散:
		for(int i=0; i<row; ++i) {
			if(fabs(x_k1[i])>=INF) {
				isDiverge=true;
				break;
			}
		}
		if(isDiverge) {
			break;
		}

	} while(fabs(diff)>=eps);

	if(isDiverge) {
		isDiverge=false;
		printf("CG发散\n");
	} else {
		printf(
		    "最终结果\n"
		    "经过%d次迭代\n"
		    "x_k1 = \n"
		    ,k);
		output1D(x_k1,row);
	}



//---------------------------------------
//第二题:
#else
	int ROW=0, row=0;
	printf("输入方阵的最大行数: ROW = ");
	scanf("%d",&ROW);
	while(row<ROW) {
		++row;
		printf("row = %d\n",row);
		iniA(a,row);
		iniB(b,row);
		eps=0.000001;
#ifdef DEBUG
		printf("矩阵A:\n");
		output2D(a,row);
		printf("系数矩阵B:\n");
		output1D(b,row);
#endif
//计算初次迭代的各个初值:

		//r=b-Ax
		memset(tempAns,0,sizeof(double)*N);		//中间结果tempAns清零
		MatrixMtimes(a,x_k0,tempAns,row);
		VectorMinus(b,tempAns,r_k0,row);

		//p=r
		VectorCopy(r_k0,p_k0,row);

		//α= r^T*r/(p^T*A*p)
		memset(tempAns,0,sizeof(double)*N);		//中间结果tempAns清零
		MatrixMtimes(a,p_k0,tempAns,row);
		alpha=VectorMtimes(r_k0,r_k0,row)/VectorMtimes(p_k0,tempAns,row);

		//x'=x+αp
		memset(tempAns,0,sizeof(double)*N);		//中间结果tempAns清零
		for(int i=0; i<row; ++i) {
			tempAns[i]=p_k0[i]*alpha;
		}
		VectorPlus(x_k0,tempAns,x_k1,row);

#ifdef DEBUG
		printf(
		    "初始值:\n"
		    "r=\n"
		);
		output1D(r_k0,row);
		printf("p = \n");
		output1D(p_k0,row);
		printf("x_k0 = \n");
		output1D(x_k0,row);
		printf("x_k1 = \n");
		output1D(x_k1,row);
		printf("alpha = %lf\n",alpha);
#endif

		int k=0;
		double diff=0.0;
		bool isDiverge=false;

		do {
			k++;

			//r=b-Ax
			memset(tempAns,0,sizeof(double)*N);		//中间结果tempAns清零
			MatrixMtimes(a,x_k1,tempAns,row);
			VectorMinus(b,tempAns,r_k0,row);
#ifdef DEBUG
			printf("tempAns = \n");
			output1D(tempAns,row);
#endif

			//β=r^TAp/p^TAp
			memset(tempAns,0,sizeof(double)*N);		//中间结果tempAns清零
			MatrixMtimes(a,p_k0,tempAns,row);
			beta=-VectorMtimes(r_k0,tempAns,row)/VectorMtimes(p_k0,tempAns,row);

			//p'=r+βp
			memset(tempAns,0,sizeof(double)*N);		//中间结果tempAns清零
			for(int i=0; i<row; ++i) {
				tempAns[i]=p_k0[i]*beta;
			}
			VectorPlus(r_k0,tempAns,p_k0,row);

			//α=r^Tp/p^TAp
			memset(tempAns,0,sizeof(double)*N);		//中间结果tempAns清零
			MatrixMtimes(a,p_k0,tempAns,row);
			alpha=VectorMtimes(r_k0,p_k0,row)/VectorMtimes(p_k0,tempAns,row);

			//x'=x+αp
			memset(tempAns,0,sizeof(double)*N);		//中间结果tempAns清零
			memset(tempAns2,0,sizeof(double)*N);	//中间结果tempAns2清零
			for(int i=0; i<row; ++i) {
				tempAns[i]=p_k0[i]*alpha;
			}
			VectorPlus(x_k1,tempAns,tempAns2,row);
			VectorCopy(tempAns2,x_k1,row);

			//diff=
			memset(tempAns,0,sizeof(double)*N);		//中间结果tempAns清零
			memset(tempAns2,0,sizeof(double)*N);	//中间结果tempAns2清零
			MatrixMtimes(a,x_k1,tempAns,row);
			VectorMinus(b,tempAns,tempAns2,row);
			diff=NormInf(tempAns2,row)/NormInf(b,row);
#ifdef DEBUG
			printf("第%d次循环\n",k);
			printf("r = \n");
			output1D(r_k0,row);
			printf("p = \n");
			output1D(p_k0,row);
			printf("x_k0 = \n");
			output1D(x_k0,row);
			printf("x_k1 = \n");
			output1D(x_k1,row);
			printf("alpha = %lf\n",alpha);
			printf("beta = %lf\n",beta);
			printf("diff = %lf\n",diff);
#endif
			//当CG的结果超过既定的最大值INF时, 判定为发散:
			for(int i=0; i<row; ++i) {
				if(fabs(x_k1[i])>=INF) {
					isDiverge=true;
					break;
				}
			}
			if(isDiverge) {
				break;
			}
		} while(fabs(diff)>=eps);

		if(isDiverge) {
			isDiverge=false;
			printf("CG发散\n");
		} else {
			printf(
			    "最终结果\n"
			    "经过%d次迭代\n"
			    "x_k1 = \n"
			    ,k);
			output1D(x_k1,row);
		}

	}


#endif

	return 0;
}

//向量乘法:
//a[1*n] * b[n*1]
double VectorMtimes(double a[N], double b[N],int n) {
	double ans=0.0;
	for(int i=0; i<n; ++i) {
		ans+=a[i]*b[i];
	}
	return ans;
}

//向量加法:
//a[n] + b[n]
void VectorPlus(double a[N],double b[N], double ans[N], int n) {
	for(int i=0; i<n; ++i) {
		ans[i]=a[i]+b[i];
	}
	return ;
}

//向量减法:
//a[n] - b[n]
void VectorMinus(double a[N],double b[N],double ans[N],int n) {
	for(int i=0; i<n; ++i) {
		ans[i]=a[i]-b[i];
	}
	return ;
}

//方阵乘法:
//a[n*n] * b[n*1]
void MatrixMtimes(double a[N][N], double b[N],double ans[N],int n) {
	memset(ans,0,sizeof(double)*N);
	for(int i=0; i<n; ++i) {
		for(int j=0; j<n; ++j) {
			ans[i]+=a[i][j]*b[j];
		}
	}
	return ;
}

//向量复制:
void VectorCopy(double src[N],double dst[N],int n) {
	for(int i=0; i<n; ++i) {
		dst[i]=src[i];
	}
	return ;
}


//求向量的2范数:
double Norm2(double a[N],int n) {
	double sum=0.0;
	for(int i=0; i<n; ++i) {
		sum+=fabs(a[i]*a[i]);
	}
	return sqrt(sum);
}

//求向量的1范数:
double Norm1(double a[N],int n) {
	double sum=0.0;
	for(int i=0; i<n; ++i) {
		sum+=fabs(a[i]);
	}
	return sum;
}

//求向量的Inf范数
double NormInf(double a[N],int n) {
	double maxValue=a[0];
	for(int i=1; i<n; ++i) {
		maxValue=maxValue<a[i]?a[i]:maxValue;
	}
	return maxValue;
}

//求N阶方阵a 的Inf范数
double MatrixNormInf(double a[N][N], int n) {
	double maxValue=Norm1(a[0],n);
	for(int i=1; i<n; ++i) {
		maxValue=maxValue<Norm1(a[i],n)?Norm1(a[i],n):maxValue;
	}
	return maxValue;
}


void iniA(double a[N][N], int n) {
	for(int i=0; i<n; ++i) {
		for(int j=0; j<n; ++j) {
			a[i][j]= 1.0/(i+j+1);
		}
	}
	return ;
}

void iniB(double b[N],int n) {
	double sum=0;
	for(int i=0; i<n; ++i) {
		sum=0;
		for(int j=0; j<n; ++j) {
			sum+=1.0/(i+j+1);
		}
		b[i]=sum;
	}
	return ;
}

void input1D(double vec[N], int n) {
	for(int i=0; i<n; ++i) {
		scanf("%lf",&vec[i]);
	}
	return ;
}

void input2D(double matrix[N][N],int n) {
	for(int i=0; i<n; ++i) {
		input1D(matrix[i],n);
	}
	return ;
}

void output1D(double vec[N],int n) {
	for(int i=0; i<n; ++i) {
		printf("%+18.10lf ",vec[i]);
	}
	putchar('\n');
	return ;
}

void output2D(double matrix[N][N], int n) {
	for(int i=0; i<n; ++i) {
		output1D(matrix[i],n);
	}
	putchar('\n');
	return ;
}

实验四:

源码:

#include <iostream>
#include <cmath>

using namespace std;

#define DEBUG

void Thomas(double a[], double b[],double c[],double f[],double ans[],int n);
double function(double x);
void output(double arr[], int n);


int main(void) {
	int n=0;
	printf("输出三次样条差值的分段数N=");
	scanf("%d",&n);

	//输入各段Xi
	double x[n+1]= {0};
	printf("输入各个段点值Xi:");
	for(int i=0; i<n+1; ++i) {
		scanf("%lf",&x[i]);
	}

	//求解Hi
	double h[n+1]= {0};
	for(int i=1; i<n+1; ++i) {
		h[i]=x[i]-x[i-1];
	}

	//输入两个端点值
	double y_0=0.0,y_n=0.0;
	printf("输出两端点导数值 Y0' 与 Yn' :");
	scanf("%lf",&y_0);
	scanf("%lf",&y_n);

	//求解Yi
	double y[n+1];
	for(int i=0; i<n+1; ++i) {
		y[i]=function(x[i]);
	}

	//几个使用到的变量:
	double d[n+1]= {0};
	double c[n+1]= {0};
	double lambda[n+1]= {0};
	double miu[n+1]= {0};
	double M[n+1]= {0};
	//进入处理循环
	for(int i=1; i<n+1; ++i) {
		d[i]=6*((y[i+1]-y[i])/h[i+1]-(y[i]-y[i-1])/h[i])/(h[i]+h[i+1]);
		c[i]=2;
		lambda[i]=h[i+1]/(h[i]+h[i+1]);
		miu[i]=1-lambda[i];
	}
	//DEBUG: d,c,lambda,miu=?

	//转化为Thomas追赶法的参数:
	double Thomas_a[n-1]= {0};
	double Thomas_b[n-1]= {0};
	double Thomas_c[n-1]= {0};
	double Thomas_f[n-1]= {0};
	double ans[n-1]= {0};
	Thomas_b[0]=c[1];
	Thomas_f[0]=d[1];
	for(int i=1; i<=n-2; ++i) {
		Thomas_a[i]=miu[i+1];
		Thomas_b[i]=c[i+1];
		Thomas_c[i-1]=lambda[i];
		Thomas_f[i]=d[i+1];
	}
	//DEBUG: Thomas_a,Thomas_b,Thomas_c,Thomas_f=?
	Thomas(Thomas_a,Thomas_b,Thomas_c,Thomas_f,ans,n-1);
	//DEBUG: ans=?

	//将追赶法的结果提取到M中:
	for(int i=1; i<n; i++) {
		M[i]=ans[i-1];
	}
	//DEBUG: M=?

	//求解端点值:
	M[0] = (6 * ((y[1] - y[0]) / h[1] - y_0) / h[1] - M[1]) / 2;
	M[n] = (6 * (y_n - (y[n] - y[n - 1]) / h[n]) / h[n] - M[n - 1]) / 2;
	//DEBUG: M=?


	//准备inputX
	int inputTime=20;
	double inputX[inputTime]= {0};		//inputX[0]为空
	for(int i=1; i<inputTime+1; ++i) {
		inputX[i]=-1.05+0.1*i;
	}
	//DEBUG: inputX=?

	

	//求解最终的S_tar(Xi)
	int tar=1;
	double S_tar[inputTime+1]= {0};
	for(int i=1; i<inputTime+1; ++i,++tar) {
		S_tar[i]=M[tar-1]*pow(x[tar]-inputX[i],3)/(6*h[tar])+
		         M[tar]*pow(inputX[i]-x[tar-1],3)/(6*h[tar])+
		         (y[tar-1]-M[tar-1]/6*h[tar]*h[tar])*(x[tar]-inputX[i])/h[tar]+
		         (y[tar]-M[tar]/6*h[tar]*h[tar])*(inputX[i]-x[tar-1])/h[tar];
	}


	//求解Lagrange差值:
	double sum=0, mul=1;
	double L_tar[inputTime+1]= {0};
	for(int i=1; i<inputTime+1; ++i) {
		sum=0;
		for(int j=0; j<n+1; ++j) {
			mul=1;
			for(int k=0; k<n+1; ++k) {
				if(k!=j) {
					mul*=(inputX[i]-x[k])/(x[j]-x[k]);
				}
			}
			sum+=mul*y[j];
		}
		L_tar[i]=sum;
	}
	//DEBUG: L_tar=?


	//输出最终结果:
	printf("\n\n最终结果:\n");
	for(int i=1; i<inputTime+1; ++i) {
		printf("x_%2d = %+.2lf    f(%+.2lf) = %+17.14lf   S_%2d = %+17.14lf   L_%2d = %+17.14lf \n",
		       i,inputX[i],inputX[i],function(inputX[i]),i,S_tar[i],i,L_tar[i]);
	}

	return 0;
}

void output(double arr[], int n) {
	for(int i=0; i<n; ++i) {
		printf("%.12lf ",arr[i]);
	}
	return ;
}

double function(double x) {
	return 1/(1+25*x*x);
}

//Thomas追赶法解n阶三对角矩阵, 注意是n阶
void Thomas(double a[], double b[],double c[],double f[],double ans[],int n) {
	double beta[n-1]= {0};
	double y[n]= {0};

	beta[0]=c[0]/b[0];
	for(int i=0; i<n; ++i) {
		beta[i]=c[i]/(b[i]-a[i]*beta[i-1]);
	}

	y[0]=f[0]/b[0];
	for(int i=1; i<n; ++i) {
		y[i]=(f[i]-a[i]*y[i-1])/(b[i]-a[i]*beta[i-1]);
	}

	ans[n-1]=y[n-1];

	for(int i=n-2; i>=0; --i) {
		ans[i]=y[i]-beta[i]*ans[i+1];
	}

	return;
}

输入输出:

输出三次样条差值的分段数N=20
输入各个段点值Xi:
-1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
输出两端点导数值 Y0’ 与 Yn’ :
0.07396449704142012

0.07396449704142012

最终结果:
x_ 1 = -0.95 f(1717986919) = +0.04244031830239 S_ 1 = +0.04242204208118 L_ 1 = -39.95244903304251
x_ 2 = -0.85 f(858993460) = +0.05245901639344 S_ 2 = +0.05243128736673 L_ 2 = +3.45495779986415
x_ 3 = -0.75 f( 0) = +0.06639004149378 S_ 3 = +0.06639405338251 L_ 3 = -0.44705196070884
x_ 4 = -0.65 f(-858993459) = +0.08648648648649 S_ 4 = +0.08647363117869 L_ 4 = +0.20242261570456
x_ 5 = -0.55 f(-1717986918) = +0.11678832116788 S_ 5 = +0.11678687446503 L_ 5 = +0.08065999342165
x_ 6 = -0.45 f(-858993460) = +0.16494845360825 S_ 6 = +0.16486455736328 L_ 6 = +0.17976262990060
x_ 7 = -0.35 f(1717986918) = +0.24615384615385 S_ 7 = +0.24626815868132 L_ 7 = +0.23844593373813
x_ 8 = -0.25 f( 0) = +0.39024390243902 S_ 8 = +0.38941957183718 L_ 8 = +0.39509305368778
x_ 9 = -0.15 f(858993460) = +0.64000000000000 S_ 9 = +0.64316893858536 L_ 9 = +0.63675533591643
x_10 = -0.05 f(-1717986912) = +0.94117647058824 S_10 = +0.93886621228293 L_10 = +0.94249037974398
x_11 = +0.05 f(-1717986912) = +0.94117647058824 S_11 = +0.93886621228293 L_11 = +0.94249037974398
x_12 = +0.15 f(858993464) = +0.64000000000000 S_12 = +0.64316893858536 L_12 = +0.63675533591643
x_13 = +0.25 f( 0) = +0.39024390243902 S_13 = +0.38941957183718 L_13 = +0.39509305368778
x_14 = +0.35 f(1717986920) = +0.24615384615385 S_14 = +0.24626815868132 L_14 = +0.23844593373813
x_15 = +0.45 f(-858993460) = +0.16494845360825 S_15 = +0.16486455736328 L_15 = +0.17976262990060
x_16 = +0.55 f(-1717986918) = +0.11678832116788 S_16 = +0.11678687446503 L_16 = +0.08065999342165
x_17 = +0.65 f(-858993458) = +0.08648648648649 S_17 = +0.08647363117869 L_17 = +0.20242261570456
x_18 = +0.75 f( 0) = +0.06639004149378 S_18 = +0.06639405338251 L_18 = -0.44705196070884
x_19 = +0.85 f(858993460) = +0.05245901639344 S_19 = +0.05243128736673 L_19 = +3.45495779986420
x_20 = +0.95 f(1717986918) = +0.04244031830239 S_20 = +0.03964837344213 L_20 = -39.95244903304321

实验五:

#include <iostream>
#include <cmath>
#include <cfloat>

//#define DEBUG

using namespace std;

//order表示拟合函数幂次
const int order=3;

double function(double x);
void Gauss_Seidel_Iteration(double ma0[order+1][order+1], double b[order+1] ,
                            double ans[order+1]) ;
void MatrixMinus(double ma0[], double ma1[], double ans[],int n);
void MatrixCopy(double src[],double dst[],int n);
double NormInf(double vec[],int n);
void Thomas(double a[], double b[],double c[],double f[],double ans[],int n) ;
void Lagrange(double L_tar[],int inputTime,double inputX[], double x[], double y[]);
void spline(int n,double x[], double y[], double inputX[],double S_tar[]) ;


int main(void) {

	//对Phi进行赋值, 这里使用函数指针和Lambda表达式
	const int MAX_ORDER=4;
	double(*Phi[MAX_ORDER+1])(double x) = {0};
	Phi[0]=[](double x)->double {return pow(x,0);};
	Phi[1]=[](double x)->double {return pow(x,1);};
	Phi[2]=[](double x)->double {return pow(x,2);};
	Phi[3]=[](double x)->double {return pow(x,3);};
	Phi[4]=[](double x)->double {return pow(x,4);};

	int max_iteration=0;
	printf("输入n的最大值max_iteration = ");
	scanf("%d",&max_iteration);
	int n=1;
	while(n<max_iteration) {
		++n;
#ifdef DEBUG		
		n=20;
#endif
		printf("第%d次迭代, n=%d",n,n);

		//测试用:
//		n=4;

		//求Xi
		//测试用:
//		double x[n+1]= {0.0,0.25,0.5,0.75,1.0};
		double x[n+1]= {0};
		for(int i=0; i<n+1; ++i) {
			x[i]=-1+2.0*i/n;
		}
		//DEBUG: X=?


		//求解Yi
		//测试用:
//		double y[n+1]= {1.0,1.2840,1.6487,2.1170,2.7183};
		double y[n+1]= {0};
		for(int i=0; i<n+1; ++i) {
			y[i]=function(x[i]);
		}
		//DEBUG: Y=?

		//求解Gram矩阵
		double gram[order+1][order+1]= {0};
		double tempSum=0.0;
		for(int i=0; i<order+1; ++i) {
			for(int j=0; j<order+1; ++j) {
				tempSum=0.0;
				for(int k=0; k<n+1; ++k) {
					tempSum+=Phi[i](x[k])*Phi[j](x[k]);
				}
				gram[i][j]=tempSum;
			}
		}
		//DEBUG: gram=?

		//求解系数矩阵d:
		double d[order+1]= {0};
		for(int i=0; i<order+1; ++i) {
			tempSum=0.0;
			for(int j=0; j<n+1; ++j) {
				tempSum+=Phi[i](x[j])*y[j];
			}
			d[i]=tempSum;
		}
		//DEBUG: d=?

		//求解线性方程组:
		double a[order+1]= {0};
		Gauss_Seidel_Iteration(gram,d,a);
		//DEBUG: a=?

		//准备inputX
		double inputX[n+1]= {0};	//inputX[0]为空
		for(int i=1; i<n+1; ++i) {
			inputX[i]=(x[i]+x[i-1])/2;
		}
		//DEBUG: inputX=?
		
		//计算最小二乘拟合的结果:
		double s_star[n+1]= {0};
		for(int i=0; i<n+1; ++i) {
			for(int j=0; j<order+1; ++j) {
				s_star[i]+=Phi[j](inputX[i])*a[j];
			}
		}


//---------------------------------------------
		//计算Lagrange差值的结果:
		double L_tar[n+1]= {0};
		if(n>2) {
			Lagrange(L_tar,n,inputX,x,y);
		}
		//DEBUG: L_tar=?
//---------------------------------------------
		//计算三次样条差值的结果:
		double S_tar[n+1]= {0};
		if(n>2) {
			spline(n,x,y,inputX,S_tar);
		}
//---------------------------------------------

		//输出结果:
		printf("\n\n最终结果:\n");
		for(int i=1; i<n+1; ++i) {
			printf("Xi = %+.3lf   S_*(Xi) = %+18.14lf   S_%2d  = %+18.14lf\n"
			       "              L_%2d    = %+18.14lf   f(Xi) = %+18.14lf\n\n"
			       ,inputX[i],s_star[i],i,S_tar[i],i,L_tar[i],function(inputX[i]));
		}
	}
	return 0;
}

double function(double x) {
	return 1/(1+25*x*x);
}

//使用 Gauss_Seidel 迭代法求解线性方程组
void Gauss_Seidel_Iteration(double ma0[order+1][order+1], double b[order+1] ,
                            double ans[order+1]) {
	double x_k1[order+1]= {0};
	double x_k0[order+1]= {0};
	double x_diff[order+1]= {0};
	double eps=0.000001;
	int k=0;
	double sum=0.0;

	do {
		k++;
		for(int i=0; i<order+1; ++i) {
			sum=0.0;
			for(int j=0; j<i; ++j) {
				sum+=ma0[i][j]*x_k1[j];
			}
			for(int j=i+1; j<order+1; ++j) {
				sum+=ma0[i][j]*x_k0[j];
			}
			sum-=b[i];
			sum/=(-ma0[i][i]);
			x_k1[i]=sum;
		}
		MatrixMinus(x_k1,x_k0,x_diff,order+1);
		MatrixCopy(x_k1,x_k0,order+1);
	} while(NormInf(x_diff,order+1)>=eps);
	MatrixCopy(x_k1,ans,order+1);
}

void MatrixMinus(double ma0[], double ma1[], double ans[],int n) {
	for(int i=0; i<n; ++i) {
		ans[i]=ma0[i]-ma1[i];
	}
	return ;
}
//向量复制 
void MatrixCopy(double src[],double dst[],int n) {
	for(int i=0; i<n; ++i) {
		dst[i]=src[i];
	}
	return ;
}
//求解向量Inf范数 
double NormInf(double vec[],int n) {
	double max=DBL_MIN;
	for(int i=0; i<n; ++i) {
		max>vec[i]?max=max:max=vec[i];
	}
	return max;
}

//三次样条差值:
//n: 分段个数
//x: 各个段点值
//y: f(x)
//inputX: 输入的测试x值
//S_tar: 最终结果
void spline(int n,double x[], double y[], double inputX[],double S_tar[]) {
	//求解Hi
	double h[n+1]= {0};
	for(int i=1; i<n+1; ++i) {
		h[i]=x[i]-x[i-1];
	}

	//输入两个端点值
	double y_0=0.07396449704142012,y_n=0.07396449704142012;
//	printf("输出两端点值 Y0' 与 Yn' :");
//	scanf("%lf",&y_0);
//	scanf("%lf",&y_n);

	//几个使用到的变量:
	double d[n+1]= {0};
	double c[n+1]= {0};
	double lambda[n+1]= {0};
	double miu[n+1]= {0};
	double M[n+1]= {0};
	//进入处理循环
	for(int i=1; i<n+1; ++i) {
		d[i]=6*((y[i+1]-y[i])/h[i+1]-(y[i]-y[i-1])/h[i])/(h[i]+h[i+1]);
		c[i]=2;
		lambda[i]=h[i+1]/(h[i]+h[i+1]);
		miu[i]=1-lambda[i];
	}
	//DEBUG: d,c,lambda,miu=?

	//转化为Thomas追赶法的参数:
	double Thomas_a[n-1]= {0};
	double Thomas_b[n-1]= {0};
	double Thomas_c[n-1]= {0};
	double Thomas_f[n-1]= {0};
	double ans[n-1]= {0};
	Thomas_b[0]=c[1];
	Thomas_f[0]=d[1];
	for(int i=1; i<=n-2; ++i) {
		Thomas_a[i]=miu[i+1];
		Thomas_b[i]=c[i+1];
		Thomas_c[i-1]=lambda[i];
		Thomas_f[i]=d[i+1];
	}
	//DEBUG: Thomas_a,Thomas_b,Thomas_c,Thomas_f=?
	Thomas(Thomas_a,Thomas_b,Thomas_c,Thomas_f,ans,n-1);
	//DEBUG: ans=?

	//将追赶法的结果提取到M中:
	for(int i=1; i<n; i++) {
		M[i]=ans[i-1];
	}
	//DEBUG: M=?

	//求解端点值:
	M[0]=(6*((y[1]-y[0])/h[1]-y_0)/h[1]-M[1])/2;
	M[n]=(6*(y_n-(y[n]-y[n-1])/h[n])/h[n]-M[n-1])/2;
	//DEBUG: M=?

	//求解最终的S20(Xi)
	int tar=1;
//	double S_tar[inputTime+1]= {0};

	for(int i=1; i<n+1; ++i,++tar) {
		S_tar[i]=M[tar-1]*pow(x[tar]-inputX[i],3)/(6*h[tar])+
		         M[tar]*pow(inputX[i]-x[tar-1],3)/(6*h[tar])+
		         (y[tar-1]-M[tar-1]/6*h[tar]*h[tar])*(x[tar]-inputX[i])/h[tar]+
		         (y[tar]-M[tar]/6*h[tar]*h[tar])*(inputX[i]-x[tar-1])/h[tar];
	}
	
	return ;
}

//---------------------------------------------
//tar: 插值多项式的次数
//L_tar: 最终结果
//inputTime: 输入的测试X的个数
//inputX[]: 输入的x的值,[inputTime+1]
//x[]: 段点值
//y[]: 段点函数值 
void Lagrange(double L_tar[],int inputTime,double inputX[], double x[], double y[]) {
	double sum=0, mul=1;
	for(int i=1; i<inputTime+1; ++i) {
		sum=0;
		for(int j=0; j<inputTime+1; ++j) {
			mul=1;
			for(int k=0; k<inputTime+1; ++k) {
				if(k!=j) {
					mul*=(inputX[i]-x[k])/(x[j]-x[k]);
				}
			}
			sum+=mul*y[j];
		}
		L_tar[i]=sum;
	}

	return ;
}

//Thomas追赶法解n阶三对角矩阵, 注意是n阶
void Thomas(double a[], double b[],double c[],double f[],double ans[],int n) {
	double beta[n-1]= {0};
	double y[n]= {0};

	beta[0]=c[0]/b[0];
	for(int i=0; i<n; ++i) {
		beta[i]=c[i]/(b[i]-a[i]*beta[i-1]);
	}

	y[0]=f[0]/b[0];
	for(int i=1; i<n; ++i) {
		y[i]=(f[i]-a[i]*y[i-1])/(b[i]-a[i]*beta[i-1]);
	}

	ans[n-1]=y[n-1];

	for(int i=n-2; i>=0; --i) {
		ans[i]=y[i]-beta[i]*ans[i+1];
	}

	return;
}

实验六:

在这里插入图片描述

源码:

#include <iostream>
#include <cmath>

//#define DEBUG

using namespace std;

//数组默认的最大容量 
const int MAX=105;

double function(double x) ;

int main(int argc, char* argv[]) {

	double a=0.0, b=0.0;
	printf("输入区间[a,b] = ");
	scanf("%lf",&a);
	scanf("%lf",&b);

	double eps=0.0;
	printf("输入精度 eps = ");
	scanf("%lf",&eps);

	double T[MAX]= {0};
	double S[MAX]= {0};
	double C[MAX]= {0};
	double R[MAX]= {0};
	double h=0.0;		//中间变量 
	double sum=0.0;		//中间变量 

	T[0]=(b-a)/2*(function(a)+function(b));

	//求到T16
	for(int temp=1; temp<5; ++temp) {
		sum=0.0;
		h=(b-a)/pow(2,temp);
		for(int i=1; i<=pow(2,temp); ++i) {
			sum+=2*function(a+h*i);
		}
		sum+=(function(a)-function(b));
		T[temp]=h/2*sum;
	}
	//求到S8
	for(int temp=0; temp<4; ++temp) {
		S[temp]=T[temp+1]+(T[temp+1]-T[temp])/3;
	}
	//求到C4
	for(int temp=0; temp<3; ++temp) {
		C[temp]=S[temp+1]+(S[temp+1]-S[temp])/15;
	}
	//求到R2:
	for(int temp=0; temp<2; ++temp) {
		R[temp]=C[temp+1]+(C[temp+1]-C[temp])/63;
	}

	//用于指向各个数组的下标
	//T[4]=T16, S[3]=S8, C[2]=C4, R[1]=R2
	int pT=4, pS=3, pC=2, pR=1;

	//最少执行2次循环 : 
	int count=0;	//用于记录循环次数 
	while(count<2||fabs(R[pR]-R[pR-1])>=eps) {
#ifdef DEBUG
		printf("%d\n",count);
#endif		
		++count;
		sum=0.0;
		h=(b-a)/pow(2,pT+1);
		for(int i=1; i<=pow(2,pT+1); ++i) {
			sum+=2*function(a+h*i);
		}
		sum+=(function(a)-function(b));
		T[pT+1]=h/2*sum;
		++pT;
		
		S[pS+1]=T[pT]+(T[pT]-T[pT-1])/3;
		++pS;
		C[pC+1]=S[pS]+(S[pS]-S[pS-1])/15;
		++pC;
		R[pR+1]=C[pC]+(C[pC]-C[pC-1])/63;
		++pR;
	}
	printf("\nT: \n");
	for(int j=0; j<=pT; ++j) {
		printf("%18.14lf", T[j]);
	}
	printf("\nS: \n");
	for(int j=0; j<=pS; ++j) {
		printf("%18.14lf", S[j]);
	}
	printf("\nC: \n");
	for(int j=0; j<=pC; ++j) {
		printf("%18.14lf", C[j]);

	}
	printf("\nR: \n");
	for(int j=0; j<=pR; ++j) {
		printf("%18.14lf", R[j]);

	}

	return 0;
}

double function(double x) {
	//特殊值,无法直接计算 
	if(!x){
		return 1;
	}
//测试用:
//	return exp(x);
	return sin(x)/x;
}

输入输出:

输入区间[a,b] = 0 1

输入精度 eps = 0.000001

T:

0.92073549240395 0.93979328480618 0.94451352166539 0.94569086358270 0.94598502993439 0.94605856096277 0.94607694306006

S:

0.94614588227359 0.94608693395179 0.94608331088847 0.94608308538495 0.94608307130556 0.94608307042583

C:

0.94608300406367 0.94608306935092 0.94608307035138 0.94608307036694 0.94608307036718

R:

0.94608307038722 0.94608307036726 0.94608307036718 0.94608307036718

Janus_V
关注
  • 1
    点赞
  • 踩
  • 2
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
专栏目录
工程计算源代码
04-11
这款软件里面的计算器实现，比较灵活，里面包含了一些工程计算，对与城建学院的学生来说，是一个很实用的工具。可以方便快捷的计算相关的公式和一些查表参数。代码注释比较详细，便于阅读，里面含有在编程时一些小技巧，你值得拥有。
工程实践代码
晴天下的大头娃娃
11-28 294
PlayerMovement using System.Collections; using System.Collections.Generic; using UnityEngine; using Cinemachine; using UnityEngine.UI; using System.Net.Sockets; using System.Net; using System.Threading; using System.Text; using System.Runtime.InteropServ
工程代码实践简单总结
一群专业码农的笔记本
08-28 1282
工程代码实践简单总结
pythonflask项目实训代码_「Python实用技巧」7.Flask工程代码架构
weixin_39546501的博客
12-16 130
【这是"冲谈技术"的第7篇原创文章】大家好，我是冲冲，一名主要从事Python服务端开发的程序猿。今天和大家分享的实用技术是“Flask工程代码架构”。首先，我们话不多说，直接上结果图：7目录结构.png看不懂？猜不透？下面开始解释：-- [目录] 公司内部的公共代码package-- [目录] 项目名称__init__.py [文件] 有此文件，整个project才是一个...
计算整个工程代码量的算法源码
qyhdt的专栏
05-24 468
去年在写软件著作权时候需要统计整体工程的代码量（由于一个大系统要写六个软件著作权），无奈之下写了个统计代码量的算法（大概写了20分钟），个人感觉很准，现在公布一下源码，基本算法就是递归也没啥高技术含量，需要用到的朋友直接拷贝 方法调用在main里面都有说明，转载时候请注明出处。 [code="java"]package com.rose.common.frame.utils; import...
计算方法课程实验报告及代码
最新发布
01-07
在本压缩包中，我们找到了一份来自合肥工业大学的计算方法课程的实验报告，以及与之相关的源代码。计算方法是计算机科学与工程、数学、物理等领域的基础课程，它研究如何利用数值方法来解决各种计算问题。以下是针对...
东北大学计算机学院软件工程实验
11-28
在本实验中，我们聚焦于东北大学计算机学院的软件工程课程，这是一门旨在培养学生软件开发、管理和维护能力的核心课程。软件工程实验是该课程的重要组成部分，通过实践性学习，学生能够深入理解软件生命周期的各个...
高性能计算：实验作业及代码
12-20
本资源“高性能计算：实验作业及代码”显然是一个与HPC相关的学习资料，包含了一些实验作业和源代码，可能是为了帮助学生或者研究人员理解和实践HPC技术。 在HPC领域，主要涉及以下几个关键知识点： 1. **并行计算...
头哥实践平台计算机组成原理计算机数据表示实验(HUST)通关代码
06-15
7. **实验代码**：实验代码通常包含对上述概念的实现，例如，编写程序来模拟二进制转换、位运算、数值表示等。通过实际操作，可以加深对这些概念的理解。 8. **满分通关**：在头哥平台上，实验可能有特定的测试用例...
数值分析（计算方法）实验代码
12-06
这个名为“数值分析（计算方法）实验代码”的压缩包文件很可能包含了多个实验项目，每个项目都对应一个具体的数值计算方法。 一、基础概念 数值分析是数学的一个分支，研究如何用计算机来近似求解数学问题。这包括...
关于机械的C程序
06-21
这是一个很齐全的计算机械齿轮的一个c 程序，比较不错哦，亲
计算方法所有实验代码+重点总结
weixin_44795952的博客
02-28 333
https://github.com/rave1sking/Computational-Methods
计算工程的代码行数
THREEK的专栏
12-18 516
看到一坨坨的代码，实在不能继续下去，有个想法，统计一下这个工程的代码量，估计一下自己的工作量。最后，很失望。 using System; using System.Collections.Generic; using System.ComponentModel; using System.Data; using System.Drawing; using System.Linq; usi
输出后的空格
qq_45733269的博客
10-21 314
int main() { int x=12; double a=3.1415926; printf("%6d,%-6d,%6d\n",x,x,x); printf("%14.10lf\n",a ); return 0; }在这里插入代码片 输出结果为 12,12 , 12 3.1415926000 %-6d整数6位但是负号左对齐 所以第二个12后还有4格 ...
《C语言及程序设计》实践项目——数据的输出
迂者-贺利坚的专栏
01-10 2147
返回：贺老师课程教学链接  C语言及程序设计初步  【阅读程序题】1、写出以下程序的输出结果，再在计算机上运行程序。对比两结果是否相同，以此检查自己的学习效果。2、遇到在视频中未讲的格式控制符，上网搜索发现其用途。（1）#include int main() { int a=12345; float b=-198.345, c=6.5; printf("a=%4d,b=%
第16节 数据的输出
m0_51439429的博客
03-11 577
一.printf函数输出 1.printf函数功能:把信息送到显示器 #include <stdio.h> int main() { printf("Hello World!\n"); printf("Welcome\nto\nBeijing!\n"); } 二.使用格式描述串 用以%开始的转换描述 转换描述 对应参数类型 %d int %ld long %f double(不是float) %lf long double %c cahr #includ
（%f）单精度与（%lf）双精度是什么意思，有什么区别？
代码教学,蓝屏警告!
07-20 9276
1.float与doubleC语言中常用的小数有两种类型，分别是float或double；float称为单精度浮点型，double称为双精度浮点型。不像整数，小数没有那么多幺蛾子，小数的长度是固定的，float始终占用4个字节，double始终占用8个字节。专业术语大家可能一时无法理解，下文采用的形式来说明1.float与double。......
%.lf（四舍五入）
热门推荐
常胜将军的博客
10-22 1万+
C语言程序中浮点数类型（%.2lf）编译器默认四舍五入，如果不需要四舍五入，则要自行处理（编写函数）。（已经写过测试代码进行了验证！！！）。
认知实习实验报告
EloflyS的博客
09-09 1092
认知实习实验报告 实验准备：实验环境的搭建 准备VMWare虚拟机平台，在官方网站上下载Kali Linux并安装在VMWare中，作为实验的攻击机；再下载一个没有修复MS17-010漏洞（永恒之蓝）的Windows7安装在虚拟机当中，并在安装完成后关闭防火墙，打开445端口，作为被攻击的靶机，到此，实验环境搭建完毕 实验一：社会工程攻击的识别与防御 实验二：信息收集和漏洞扫描 为了对靶机进行攻击，我们首先需要获取该靶机的IP地址，获取方式就是在靶机中打开终端cmd，在其中输入以下命令 ipconfig
MATLAB控制工程实验：模拟典型环节
实验过程中，学生需要对每个环节的MATLAB代码进行注释，并撰写实验报告，以深入理解仿真结果和系统特性。这样的实验有助于培养学生的理论知识与实践技能，使他们能够在实际问题中应用控制理论。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值