project.cpp

这是数值分析的作业

#include <stdio.h>
#include <math.h>
#include <iostream.h>
#include <iomanip.h>
#include <stdlib.h>
#include <fstream.h>


int min(int x,int y){
	return (x >= y)?y:x;
}

int max(int x,int y){
	return (x >= y)?x:y;
}
//生成对称矩阵的压缩存储数组a[5][501]
void Create_A(double a[5][501]){
	for (int i = 1;i <= 501;i++) {
		a[0][i-1] = a[4][i-1] = -0.064;
		a[1][i-1] = a[3][i-1] = 0.16;
		a[2][i-1] = (1.64 - 0.024 * i) * sin(0.2 * i) - 0.64 * exp(0.1 / i);
	}
	a[0][0] = a[0][1] = a[1][0] = a[3][500] = a[4][499] = a[4][500] = 0;
}

//LU分解
void LUAnalysis(double a[5][501])
{
	int i,j,t;
	for(i = 0; i < 501; i++)
	{
		for(j = i; j <= min(i + 2,500); j++)
		{
			for(t = max(0,max(i - 2, j - 2)); t <= i - 1; t++)
				a[i - j + 2][j] -= a[i - t + 2][t] * a[t - j + 2][j];
		}
		for(j = i + 1; j <= min(i + 2,500); j++)
		{
			for(t = max(0,max(i - 2, j - 2)); t <= i - 1; t++)
				a[j - i + 2][i] -= a[j - t + 2][t] * a[t - i + 2][i];
			a[j - i + 2][i] /= a[2][i];
		}
	}
}

//LU分解后求解方程组 (LU)x = b
void LUSolution(double a[5][501],double b[501],double x[501])
{
	int i,t;
	double temp;
	for(i = 0; i < 501; i++)
	{
		x[i] = b[i];
		for(t = max(0,i-2); t < i; t++)
			x[i] -= a[i-t+2][t] * x[t];
	}
	for(i = 500; i >= 0; i--)
	{
		temp = 0;
		for(t = i+1; t <= min(i+2,500); t++)
			temp += x[t] * a[i-t+2][t];
		x[i] -= temp;
		x[i] /= a[2][i];
	}
}


//用幂法求最大特征值
double Power(double a[5][501],double OFFSET) {
	int temp_min,i,j;
	Create_A(a);
	for (i=0;i<501;i++) {
		a[2][i] -= OFFSET;
	}
	double u[501] = {0};
	for (int q=0;q<501;q++) {
		u[q]=1;
	}
	double y[501] = {0};
	double delta,sum,nn,bb1,bb2;
	bool first = true;
	delta = 1;
	while (delta > 1e-12) {
		sum = 0;
		for (i=1;i<=501;i++) {
			sum +=u[i-1]*u[i-1];
		}
		nn = sqrt(sum);
		for (i=1;i<=501;i++) {
			y[i-1] = u[i-1]/nn;
		}
		for (i=1;i<=501;i++) {
			u[i-1] = 0;
			temp_min = min(501,i+2);
			for (j=max(1,i-2);j<=temp_min;j++) {
				u[i-1] = u[i-1] + a[i-j+2][j-1]*y[j-1];
			}
		}
		if (first == true) {
			bb1 = 0;
			for (i=1;i<=501;i++) {
				bb1 += y[i-1]*u[i-1];
			}
			first = false;
		}
		else{
			bb2 = 0;
			for (i=1;i<=501;i++) {
				bb2 += y[i-1]*u[i-1];
			}
			delta = fabs(bb2-bb1)/fabs(bb2);
			bb1 = bb2;
		}
	}
	return bb2;
}

//用反幂法求最小特征值
double AntiPower(double a[5][501],double OFFSET){
	int i;
	Create_A(a);
	for (i=0;i<501;i++) {
		a[2][i] -= OFFSET;
	}
	double delta,sum,nn,bb1,bb2;
	bool flag = true;
	delta = 1;
	double u[501] = {0};
	for (i=0;i<501;i++) {
		u[i]=1.0;
	}
	double y[501] = {0};
	LUAnalysis(a);
	while (delta > 1e-12) {
		sum = 0;
		for (i=0;i<501;i++) {
			sum += u[i]*u[i];
		}
		nn = sqrt(sum);
		for (i=0;i<501;i++) {
			y[i] = u[i]/nn;
		}
		LUSolution(a,y,u);
		if (flag == true) {
			bb1 = 0;
			for (i=0;i<501;i++) {
				bb1 += y[i]*u[i];
			}
			flag = false;
		}
		else{
			bb2 = 0;
			for (i=0;i<501;i++) {
				bb2 += y[i]*u[i];
			}
			delta = fabs(1/bb2-1/bb1)/fabs(1/bb2);
			bb1 = bb2;
		}
	}
	return 1/bb2;
}

void main(){
	double A[5][501],AA[5][501];
	double offset,u,det;
	double AbsMaxCharacter,AbsMinCharacter,MaxCharacter,MinCharacter,TempCharacter;
	ofstream Res("d:\\myresult.txt");
	Res<<setprecision(12)<<setiosflags(ios::scientific);
	AbsMaxCharacter = Power(A,0);
	if(AbsMaxCharacter > 0)	//按模最大特征值即是最大特征值
	{
		MaxCharacter = AbsMaxCharacter;
		Res<<"最大特征值为："
			<<MaxCharacter<<endl;;
		offset =fabs(MaxCharacter);
		MinCharacter = Power(A,offset) + offset;
		Res<<"最小特征值为："
			<<MinCharacter<<endl;
	}
	else	//按模最大特征值即是最小特征值
	{
		MinCharacter = AbsMaxCharacter;
		Res<<"最小特征值为："
			<<MinCharacter<<endl;;
		offset = MinCharacter;
		MaxCharacter = Power(A,offset) + offset;
			Res<<"最大特征值为："
			<<MaxCharacter<<endl;
	}
	AbsMinCharacter = AntiPower(A,0);
	Res<<"模最小的特征值为："
		<<AbsMinCharacter<<endl;
	for (int k=1;k<40;k++) {
		u = MinCharacter +  (MaxCharacter - MinCharacter) *(double(k) / 40.0);
		Res<<"与u"<<k<<"（"<<u<<"）最接近的特征值为："
			<<AntiPower(A,u) + u<<endl;
	}
	Res<<"条件数cond(A)2为："
		<<fabs(AbsMaxCharacter/AbsMinCharacter)<<endl;
	Create_A(A);
	LUAnalysis(A);
	det = 1;
	for (int i=0;i<501;i++) {
		det *= A[2][i];
	}
	Res<<"行列式detA为："
		<<det;
}