#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;
}
C++DEBUG
最新推荐文章于 2024-07-21 13:36:11 发布