Numerical Analysis Computer Assignment First Question
#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;
}
//Initialize U to be a vector of all ones
for (int i = 0; i < n; i + + ) {<!-- -->
U[i] = 1.0;
// if (i==500){<!-- -->
// U[i] = 1.0;
// }
}
while (enow >= e) {<!-- -->
//Normalized vector 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\\
", y[i]);
}
// matrix-vector multiplication 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];
}
}
// Calculate new beta
double betanow = 0.0;
for (int i = 0; i < n; i + + ) {<!-- -->
betanow + = y[i] * U[i];
}
// Calculate enow
enow = fabs((betanow - beta) / betanow);
beta = betanow;
TIMES + = 1;
}
//printf("This iteration used step: %d\\
", TIMES);
//printf("%d\\
", 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 + + ) {<!-- -->
// The first for loop iterates each column, i represents the current column
for (int k = i; k < 501; k + + ) {<!-- -->
//The elements of U are equal to the corresponding elements of A (upper triangle part), copy directly
U[i][k] = A[i][k];
}
//The second for loop iterates each row, k represents the current row
for (int k = i; k < 501; k + + ) {<!-- -->
//The elements of L are equal to the corresponding elements of A divided by the diagonal elements of U
L[k][i] = A[k][i] / U[i][i];
}
//The third for loop iterates the following rows
for (int j = i + 1; j < 501; j + + ) {<!-- -->
//Iterate through each column (k represents the current column)
for (int k = i + 1; k < 501; k + + ) {<!-- -->
// The off-diagonal elements of A are equal to the original value of A minus the corresponding element of L multiplied by the corresponding element of U
A[j][k] = A[j][k] - L[j][i] * U[i][k];
}
}
}
//Initialize x to a vector of all ones
for (int i = 0; i < n; i + + ) {<!-- -->
x[i] = 1.0;
// if (i==500){<!-- -->
// x[i] = 1.0;
// }
}
while (enow >= e) {<!-- -->
// normalized vector 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;
}
//Solve 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];
}
// Calculate new beta
double betanow = 0.0;
for (int i = 0; i < n; i + + ) {<!-- -->
betanow + = y[i] * x[i];
}
// Calculate enow
enow = fabs((betanow - beta) / betanow);
beta = betanow;
TIMES + = 1;
}
//printf("This iteration used step: %d\\
", TIMES);
//printf("%d\\
", 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 + + ) {<!-- -->
// The first for loop iterates each column, i represents the current column
for (int k = i; k < 501; k + + ) {<!-- -->
//The elements of U are equal to the corresponding elements of A (upper triangle part), copy directly
U[i][k] = A[i][k];
}
//The second for loop iterates each row, k represents the current row
for (int k = i; k < 501; k + + ) {<!-- -->
//The elements of L are equal to the corresponding elements of A divided by the diagonal elements of U
L[k][i] = A[k][i] / U[i][i];
}
//The third for loop iterates the following rows
for (int j = i + 1; j < 501; j + + ) {<!-- -->
//Iterate through each column (k represents the current column)
for (int k = i + 1; k < 501; k + + ) {<!-- -->
// The off-diagonal elements of A are equal to the original value of A minus the corresponding element of L multiplied by the corresponding element of 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;
//Define A matrix
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\\
", lamda1);
lamda501 = mifa(A,lamda1);
printf("lamda502 is: %.12e\\
", lamda501);
for(int i = 0;i<39;i + + ){<!-- -->
uk[i + 1] = lamda1 + (i + 1)*(lamda501-lamda1)/40;
}
for ( int i = 0;i<40;i + + ){<!-- -->
//Define A matrix
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\\
", lamda);
}
cond2 = fabs(lamda1)/fabs(lamdas);
printf("This condiction number2 is: %.12lf\\
", cond2);
//Define A matrix
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\\
", detA);
return 0;
}```