数值分析-幂法和反幂法C语言

数值分析上机作业第一题

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

double mifa(double A[501][501],double pianyi){
	int i,rows = 501,cols = 501,n = 501,TIMES=0;
	double U[501],y[501],e = 1e-12,enow = DBL_MAX,beta = DBL_MAX;
	for (int i = 0; i < n; i++) {
        A[i][i] = A[i][i]-pianyi;
    }
	//初始化 U 为一个全为1的向量
    for (int i = 0; i < n; i++) {
        U[i] = 1.0;
//        if (i==500){
//    		U[i] = 1.0;
//		}
    }
    
    while (enow >= e) {
        // 归一化向量 U
        double norm_U = 0.0;
        for (int i = 0; i < n; i++) {
            norm_U += U[i] * U[i];
        }
        norm_U = sqrt(norm_U);
        for (int i = 0; i < n; i++) {
            y[i] = U[i] / norm_U;
            //printf("beta: %lf\n", y[i]);
        }
        
        // 矩阵-向量相乘 A*y
        for (int i = 0; i < n; i++) {
            U[i] = 0.0;
            for (int j = 0; j < n; j++) {
                U[i] += A[i][j] * y[j];
            }
        }

        // 计算新的 beta
        double betanow = 0.0;
        for (int i = 0; i < n; i++) {
            betanow += y[i] * U[i];
        }

        // 计算 enow
        enow = fabs((betanow - beta) / betanow);
        beta = betanow;
        
        TIMES += 1;
    }
    //printf("This iteration used step: %d\n", TIMES);
    //printf("%d\n", TIMES);
    return beta+pianyi;
}

double fanmifa(double A[501][501],double pianyi) {
	double L[501][501]={0},U[501][501]={0},norm_x,x[501],n=501,tempy[501]={0},e = 1e-12,enow = DBL_MAX,beta = DBL_MAX,y[501];
	int TIMES=0;
	for (int i = 0; i < n; i++) {
        A[i][i] = A[i][i]-pianyi;
    }
    for (int i = 0; i < 501; i++) {
        // 第一个for循环迭代每一列，i表示当前列
        for (int k = i; k < 501; k++) {
            // U的元素等于A的对应元素（上三角部分），直接复制
            U[i][k] = A[i][k];
        }
        
        // 第二个for循环迭代每一行，k表示当前行
        for (int k = i; k < 501; k++) {
            // L的元素等于A的对应元素除以U的对角元素
            L[k][i] = A[k][i] / U[i][i];
        }
        
        // 第三个for循环迭代下面的行
        for (int j = i + 1; j < 501; j++) {
            // 迭代每一列（k表示当前列）
            for (int k = i + 1; k < 501; k++) {
                // A的非对角元素等于A的原值减去L的对应元素乘以U的对应元素
                A[j][k] = A[j][k] - L[j][i] * U[i][k];
            }
        }
    }
	//初始化 x 为一个全为1的向量
    for (int i = 0; i < n; i++) {
    	x[i] = 1.0;
//        if (i==500){
//    		x[i] = 1.0;
//		}
    }
    while (enow >= e) {
        // 归一化向量 x
        double norm_x = 0.0;
        for (int i = 0; i < n; i++) {
            norm_x += x[i] * x[i];
        }
        norm_x = sqrt(norm_x);
        for (int i = 0; i < n; i++) {
            y[i] = x[i] / norm_x;
        }
        
        // 求解uk 
        for (int i = 0; i < n; i++) {
        	tempy[i] = y[i];
        	for (int j = 0; j < i; j++) {
            	tempy[i] -= L[i][j] * tempy[j];
        	}
    	}
		for (int i = n - 1; i >= 0; i--) {
        	x[i] = tempy[i];
        	for (int j = i + 1; j < n; j++) {
            	x[i] -= U[i][j] * x[j];
        	}
        	x[i] /= U[i][i];
    	}

        // 计算新的 beta
        double betanow = 0.0;
        for (int i = 0; i < n; i++) {
            betanow += y[i] * x[i];
        }

        // 计算 enow
        enow = fabs((betanow - beta) / betanow);
        beta = betanow;
        TIMES += 1;
    }
    //printf("This iteration used step: %d\n", TIMES);
    //printf("%d\n", TIMES);
    return 1/(beta)+pianyi;
}

double det(double A[501][501]) {
	double L[501][501]={0},U[501][501]={0},detA = 1;
    for (int i = 0; i < 501; i++) {
        // 第一个for循环迭代每一列，i表示当前列
        for (int k = i; k < 501; k++) {
            // U的元素等于A的对应元素（上三角部分），直接复制
            U[i][k] = A[i][k];
        }
        
        // 第二个for循环迭代每一行，k表示当前行
        for (int k = i; k < 501; k++) {
            // L的元素等于A的对应元素除以U的对角元素
            L[k][i] = A[k][i] / U[i][i];
        }
        
        // 第三个for循环迭代下面的行
        for (int j = i + 1; j < 501; j++) {
            // 迭代每一列（k表示当前列）
            for (int k = i + 1; k < 501; k++) {
                // A的非对角元素等于A的原值减去L的对应元素乘以U的对应元素
                A[j][k] = A[j][k] - L[j][i] * U[i][k];
            }
        }
    }
    for (int i = 0; i < 501; i++) {
    	detA *= U[i][i];
    }
    return detA;
}
int main() {
    double A[501][501] = {0};
    double e = 1e-12,lamda1,lamda501,lamdas,uk[40],cond2,lamda,detA;
    int ITERS;
    int n = 501;

    //定义A矩阵 
    for (ITERS = 1; ITERS < 502; ITERS++) {
        A[ITERS-1][ITERS-1] = (1.64 - 0.024 * ITERS) * sin(0.2 * ITERS) - 0.64 * exp(0.1 / ITERS);
        if (ITERS > 1) {
            A[ITERS - 2][ITERS-1] = 0.16;
            A[ITERS-1][ITERS - 2] = 0.16;
        }
        if (ITERS > 2) {
            A[ITERS - 3][ITERS-1] = -0.064;
            A[ITERS-1][ITERS - 3] = -0.064;
        }
    }
    
    lamda1 = mifa(A,0);
    printf("lamda1 is: %.12e\n", lamda1);
    lamda501 = mifa(A,lamda1);
    printf("lamda502 is: %.12e\n", lamda501);
    for(int i = 0;i<39;i++){
    	uk[i+1] = lamda1+(i+1)*(lamda501-lamda1)/40;
	}
	for ( int i = 0;i<40;i++){
		//定义A矩阵 
    	for (ITERS = 1; ITERS < 502; ITERS++) {
        	A[ITERS-1][ITERS-1] = (1.64 - 0.024 * ITERS) * sin(0.2 * ITERS) - 0.64 * exp(0.1 / ITERS);
        	if (ITERS > 1) {
            	A[ITERS - 2][ITERS-1] = 0.16;
            	A[ITERS-1][ITERS - 2] = 0.16;
        	}
        	if (ITERS > 2) {
            	A[ITERS - 3][ITERS-1] = -0.064;
            	A[ITERS-1][ITERS - 3] = -0.064;
        	}
    	}
    	lamda = fanmifa(A,uk[i]);
    	if (i == 0){
    		lamdas = lamda;
		}
    	printf("This step lamda: %.12e\n", lamda);
	}
    
    cond2 = fabs(lamda1)/fabs(lamdas);
    printf("This condiction number2 is: %.12lf\n", cond2);
    //定义A矩阵 
    	for (ITERS = 1; ITERS < 502; ITERS++) {
        	A[ITERS-1][ITERS-1] = (1.64 - 0.024 * ITERS) * sin(0.2 * ITERS) - 0.64 * exp(0.1 / ITERS);
        	if (ITERS > 1) {
            	A[ITERS - 2][ITERS-1] = 0.16;
            	A[ITERS-1][ITERS - 2] = 0.16;
        	}
        	if (ITERS > 2) {
            	A[ITERS - 3][ITERS-1] = -0.064;
            	A[ITERS-1][ITERS - 3] = -0.064;
        	}
    	}
    detA = det(A);
    printf("det of A is: %.12e\n", detA);
    return 0;
}```