C++DEBUG

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<iostream>
#include<vector>
using namespace std;

/定义函数/

//最大最小值函数 返回数temp
int min(int a, int b)
{
	int temp;
	if (a >= b)
	{
		temp = b;
	}
	else
	{
		temp = a;
	}
	return temp;
}
int max(int a, int b)
{
	int temp;
	if (a >= b)
	{
		temp = a;
	}
	else
	{
		temp = b;
	}
	return temp;
}

//矩阵、向量运算函数/

///二范数、向量的模 返回数ans
double vector_mo(double k[])
{
	double sum = 0;
	for (int i = 0; i < 501; i++)
	{
		sum += k[i] * k[i];
	}
	double ans = sqrt(sum);
	return ans;
}

///向量内积 返回数sum
double vector_vector(double a[], double b[])//两向量乘，默认a、b元素个数相等
{
	double sum = 0;
	for (int i = 0; i < sizeof(a); i++)
	{
		sum += a[i] * b[i];
	}
	return sum;
}

//得到幂法中的beta
double get_beta(double y[], double u[])
{
	double B = 0;
	for (int i = 0; i < 501; i++)
	{
		B += y[i] * u[i];
	}
	return B;
}

//A矩阵右乘向量
double *matrix_vector(vector<vector<double>> matrixD, double vectorD1[501])  //不是用A矩阵而是C矩阵进行运算
{


	double *A;
	A = (double*)malloc(501 * sizeof(double));

	A[0] = matrixD[2][0] * vectorD1[0] + matrixD[1][1] * vectorD1[1] + matrixD[0][2] * vectorD1[2];
	A[1] = matrixD[3][0] * vectorD1[0] + matrixD[2][1] * vectorD1[1] + matrixD[1][2] * vectorD1[2] + matrixD[0][3] * vectorD1[3];
	A[499] = matrixD[4][497] * vectorD1[497] + matrixD[3][498] * vectorD1[498] + matrixD[2][499] * vectorD1[499] + matrixD[1][500] * vectorD1[500];
	A[500] = matrixD[4][498] * vectorD1[498] + matrixD[3][499] * vectorD1[499] + matrixD[2][500] * vectorD1[500];
	for (int i = 2; i < 499; i++)
	{
		double sum2 = 0;
		for (int j = -2; j <= 2; j++)
		{
			sum2 += matrixD[2 + j][i] * vectorD1[i + j];
		}
		A[i] = sum2;
	}
	return A;//得到uk向量
}


///A矩阵的逆右乘向量///
double *anti_matrix_vector(double matrix[5][501], double vector[501])
{
	double sum3, sum4;
	double *x;  //存储uk
	x = (double*)malloc(501 * sizeof(double));
	double e[501];  //存储被赋值后的向量，不改变原向量

	e[0] = vector[0];
	for (int i = 2; i <= 501; i++)
	{
		sum3 = 0;
		for (int t = max(1, i - 2); t <= i - 1; t++)
		{
			sum3 = sum3 + matrix[i - t + 2][t - 1] * e[t - 1];
		}
		e[i - 1] = vector[i - 1] - sum3;
	}

	x[500] = e[500] / matrix[2][500];

	for (int i = 500; i >= 1; i--)
	{
		sum4 = 0;
		for (int t = i + 1; t <= min(i + 2, 501); t++)
		{
			sum4 = sum4 + matrix[i - t + 2][t - 1] * x[t - 1];
		}

		x[i - 1] = (e[i - 1] - sum4) / matrix[2][i - 1];
	}
	return x;
}


///幂法///
double mifa(vector<vector<double>> a, double u0[])
{

	double *uk;
	uk = (double*)malloc(501 * sizeof(double));//用于存储迭代得到的uk
	double *y;
	y = (double*)malloc(501 * sizeof(double));//用于存储迭代得到的uk
	for (int i = 0; i<501; i++)
	{
		uk[i] = u0[i];
	}                                         //不改变输入向量，用新向量代替
	double beta1;//让程序运行，并存储betak
	double beta0;//当做beta_k-1
	double temp = 0;
	do
	{
		//将上一次迭代得到的beta_k作为这一次迭代的beta_k-1
		beta0 = temp;

		//获得ita值
		double ita = vector_mo(uk);

		//得到y向量
		for (int i = 0; i < 501; i++)
		{
			y[i] = uk[i] / ita;
		}

		//得到uk向量
		uk = matrix_vector(a, y);

		//得到beta
		beta1 = get_beta(y, uk);

		temp = beta1;

	} while (fabs(beta1 - beta0) / fabs(beta1) > 1e-12);
	double lamda = beta1;

	return  lamda;

}

///反幂法///
double fanmifa(double a[][501], double u0[])
{

	double *uk1;
	uk1 = (double*)malloc(501 * sizeof(double));//用于存储迭代得到的uk
	double *y1;
	y1 = (double*)malloc(501 * sizeof(double));//用于存储迭代得到的uk

	for (int i = 0; i<501; i++)
	{
		uk1[i] = u0[i];

	}
	double beta1;//让程序运行，并存储betak
	double beta0;//当做beta_k-1
	double temp = 0;



	do
	{
		//将上一次迭代得到的betak作为这一次迭代的k-1
		beta0 = temp;

		//获得ita值
		double ita = vector_mo(uk1);

		//得到y向量
		for (int i = 0; i < 501; i++)
		{
			y1[i] = uk1[i] / ita;
		}

		//得到uk向量
		uk1 = anti_matrix_vector(a, y1);

		//得到beta
		beta1 = get_beta(y1, uk1);
		beta1 = 1 / beta1;
		temp = beta1;

	} while (fabs(beta1 - beta0) / fabs(beta1) > 1e-12);

	double lamda = beta1;

	return  lamda;

}

///LU分解，元素存储在原矩阵中///
void Doolittle_LU(double matrix[][501])
{

	double sum1, sum2;  //定义中间量，不在循环中改变原矩阵
	for (int k = 1; k <= 501; k++)
	{
		for (int j = k; j <= min(k + 2, 501); j++)
		{
			sum1 = 0;
			for (int t = max(1, max(k - 2, j - 2)); t <= (k - 1); t++)
			{
				sum1 = sum1 + matrix[k - t + 2][t - 1] * matrix[t - j + 2][j - 1];
			}
			matrix[k - j + 2][j - 1] = matrix[k - j + 2][j - 1] - sum1;

		}

		if (k == 501)
		{
			break;
		}

		for (int i = k + 1; i <= min(k + 2, 501); i++)
		{
			sum2 = 0;
			for (int t = max(1, max(k - 2, i - 2)); t <= (k - 1); t++)
			{
				sum2 = sum2 + matrix[i - t + 2][t - 1] * matrix[t - k + 2][k - 1];
			}
			matrix[i - k + 2][k - 1] = (matrix[i - k + 2][k - 1] - sum2) / matrix[2][k - 1];
		}
	}
}



///定义函数结束/

//主程序
int main()
{
    设定初始向量u0
    double u0[501];
    double b = 0.16;
    double c = -0.064;
    double k1, k2, temp, lamda_1, lamda_501, lamda_s;//k1,k2:幂法得到的值;temp:中间量;lamda_1, lamda_501:得到的最终结果
    ///创建初始矩阵C///

    vector<vector<double>> C(5); //这里m就是相当于二维数组的行数，必须写，不然报错
    //这里创建一个m*n的二维vector
    for(int i = 0 ; i < 5; i++){//这里是给内层vector定义大小。默认是0,这里n是个数，不是值
        C[i].resize(501);  //利用resize进行扩充
    }

    ///储存矩阵A元素的数组C,用于三角分解法解带状线性方程组///
    C[0][0] = 0;
    C[0][1] = 0;
    C[1][0] = 0;
    C[3][500] = 0;
    C[4][500] = 0;
    C[4][499] = 0;//C中零元素
    for (int i = 0; i < 501; i++)//i 0-500
    {
        double num1 = 0.1 / (i + 1);
        C[2][i] = double((1.64 - 0.024*(i + 1))*sin(0.2*(i + 1))) - 0.64*exp(num1);
        if (i > 0)//储存b值
        {
            C[1][i] = b;
            C[3][i - 1] = b;
        }
        if (i > 1)//储存c值
        {
            C[0][i] = c;
            C[4][i - 2] = c;
        }
    }

    ///赋值初始u0///
    for (int i = 0; i < 501; i++)
    {
        u0[i] = 1;
    }
    u0[500] = 0;

    //使用幂法计算lamda_1和lamda_501//

    //k1,k2:幂法得到的值;temp:中间量;lamda_1, lamda_501:得到的最终结果
    //对初始矩阵使用幂法
    k1 = mifa(C, u0);
    //原点平移，得到B矩阵//
    for (int i = 0; i < 501; i++)
    {
        C[2][i] = C[2][i] - k1;
    }
    //对平移后矩阵使用幂法//
    k2 = mifa(C, u0);
    //得到limda_0//
    temp = k1 + k2;
    //判断limda_0与k1大小，大值为lamda_501，小值为lamda_1
    if (temp > k1)
    {
        lamda_501 = temp;
        lamda_1 = k1;
    }
    else
    {
        lamda_1 = temp;
        lamda_501 = k1;
    }

    ///幂法计算lamda_1和lamda_501结束，输出相关值///
    printf("lamda_1为%12.11e\n", lamda_1);
    printf("lamda_501为%12.11e\n", lamda_501);

    程序结束
    system("pause");
    return 0;
}