역 거듭제곱

언어/Coding 2015. 12. 22. 02:47

/*LU분해 후 역행렬을 구하고 난뒤 최대 고유값을 구한다*/

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

double A[10][10],x0[10],sum[10],xR[10],B[10][10],c[10][10],L[10][10],U[10][10],b2[10],E[10],x[10],B[10][10];

int n;

void LUCrout();

void Reverse();


void LUCrout() {

        int i,j,k;

        double sum;

        for(i=1;i<=n;i++)

                L[i][1]=A[i][1];

        for(j=1;j<=n;j++)

                U[1][j]=A[1][j]/L[1][1];

        for(j=2;j<=n;j++) {

                for(i=j;i<=n;i++) {

                        sum=0;

                        for(k=1;k<=j-1;k++) sum=sum+L[i][k]*U[k][j];

                        L[i][j]=A[i][j]-sum;

                }

                U[j][j]=1;

                for(i=j+1;i<=n;i++) {

                        sum=0.0;

                        for(k=1;k<=j-1;k++) sum=sum+L[j][k]*U[k][i];

                        U[j][i]=(A[j][i]-sum)/L[j][j];

                }

        }

}//end of LuCrout


void Reverse() {

        int i,j,a;

        double sum;

        printf("역행렬을 구합니다\n");

for(a=1;a<=n;a++){

        printf("역행렬을 구합니다. 단위 행렬을 입력해주세요\n");

        for(i=1;i<=n;i++){

                printf("E[%d]:",i);

                scanf("%lf", &E[i]);

        }


        for(i=1;i<=n;i++){

                sum=0;

                for(j=1;j<=i-1;j++)

                        sum=sum+L[i][j]*b2[j];

                b2[i]=(E[i]-sum)/L[i][i];

        }

        printf("\n");


        for(i=n;i>=1;i--){

                sum=0;

                for(j=i+1;j<=n;j++)

                        sum=sum+U[i][j]*x[j];

                x[i]=b2[i]-sum;

        }

        for(i=1;i<=n;i++)

                printf("역행렬 x[%d]:%lf\n",i,x[i]);

        for(i=1;i<=n;i++)

                B[i][a]=x[i];

}

}// end of Reverse


int main(){


        int i,j,k,repeat,v;

        float ramda,temp,temp2,x1,deter;

        printf("몇 행렬: ");

        scanf("%d",&n);

printf("반복 횟수 입력:");

        scanf("%d",&repeat);

        printf("기본 행렬을 입력해 주세요:\n");

        for(i=1;i<=n;i++) //행 입력

        {

                for(j=1;j<=n;j++) //열 입력

                {

                        printf("A[%d][%d]:",i,j);

                        scanf("%lf",&A[i][j]);

                }

        }

        LUCrout();

        Reverse();

        printf("역행렬(B): \n");

        for(i=1;i<=n;i++)

        {

                for(j=1;j<=n;j++)

                {

                        printf("%f\t",B[i][j]);

                }

                printf("\n");

        }

        for(i=1;i<=n;i++)

        {

                printf("초기 x(0)을 입력해 주세요: ");

                scanf("%lf",&x0[i]);

        }

        v=1;

        for(i=1;i<=repeat;i++) //수렴기준 대신에 반복 횟수로 카운트

        {

                for(j=1;j<=n;j++)

                {

                        sum[j]=0;

                        for(k=1;k<=n;k++)

                        {


sum[j]=sum[j]+B[j][k]*x0[k]; // Y = sum, A역행렬=B

                        }

                        //printf("Y값 계산:%f\n",sum[j]);

                }

                temp = 0;

                temp2 = 0;

                for(j=1;j<=n;j++)

                {

                        temp=x0[j]*sum[j]+temp; // temp는 lamda구하기 위한 임시 변수

                        //printf("x0[%d]=%f, temp:%f\t",j,x0[j], temp);

                        temp2=x0[j]*x0[j]+temp2;

                        //printf("temp2:%f\n",temp2);

                        ramda=temp/temp2;

                }

                //printf("Lamda:%f\n",ramda);

                x1=0;

                for(j=1;j<=n;j++)

                        x1=pow(sum[j],2)+x1;

                x1=sqrt(x1);

                //printf("x(1):%f\n",x1);

                for(j=1;j<=n;j++)

                        x0[j]=sum[j]/x1;

                printf("반복횟수:%d\t Lamda2:%f\t Lamda:%f\t x(1):%f\t x(2):%f\t x(3):%f\n",v,ramda,1/ramda,x0[1],x0[2],x0[3]);

                v++;

        }//end of repeat

}//end of main



'언어 > Coding' 카테고리의 다른 글

simpson 1/3 공식  (0) 2015.12.22
거듭제곱  (0) 2015.12.22
Lagrange  (0) 2015.12.22
Gauss Seidel  (0) 2015.12.21
Jacobi  (0) 2015.12.21
Posted by 知彼知己百戰不殆
,

Lagrange

언어/Coding 2015. 12. 22. 02:45

#include <stdio.h>

#include <math.h>

#define nd 5

double xd[nd]={0,1,2,4,5};

double yd[nd]={0,16,48,88,0};

double Lagrange(double x);


int main() {

        double x=3.0,y;

        y=Lagrange(x);

        printf("\n x=%f : f(x)=%f\n",x,y);

 }

double Lagrange(double x)

{

        double p, sum;

        int n,i,j;

        n=nd-1;

        sum=0;

        for(i=0;i<=n;i++) {

                p=1;

                for(j=0;j<=n;j++) {

                        if(i!=j)

                        p=p*(x-xd[j])/(xd[i]-xd[j]);

                        printf("p값:%f\n",p);

                }

                sum=sum+p*yd[i];

        }

        return(sum);

}



'언어 > Coding' 카테고리의 다른 글

거듭제곱  (0) 2015.12.22
역 거듭제곱  (0) 2015.12.22
Gauss Seidel  (0) 2015.12.21
Jacobi  (0) 2015.12.21
GaussN  (0) 2015.12.21
Posted by 知彼知己百戰不殆
,

Gauss Seidel

언어/Coding 2015. 12. 21. 01:04

#include <stdio.h>

#include <math.h>

/* 대각 지배행렬이어야 성립하므로 방정식이 대각 지배 행렬에 안 맞을 경우 바꿔줘야 함  */

#define nd 20


double x[nd],a[nd][nd], b[nd];

int imax = 1000, n, iter;

double epsilon = 0.00001, resavg;

void Jacobi();


int main(){

        int i, j;


        printf(" \n 방정식 수 = ");

        scanf_s("%d", &n);

        printf(" 계수 행령a(i,j)을 입력하시오 \n");

        for (i = 1; i <= n; i = i + 1){

                printf("a(%d,j) j=1,,,n =", i);

                for (j = 1; j <= n; j = j + 1) scanf("%lf", &a[i][j]);

        }

        printf(" 상수행렬 b[i]를 입력하시오 \n");

        printf("b[i] i=1,,,n = ");

        for (i = 1; i <= n; i = i + 1)scanf("%lf", &b[i]);

        printf(" 가정값 x0[i] 입력하시오 \n");

        printf(" x0[i] i=1,,,n = ");

        for (i=1;i<=n;i=i+1)scanf("%lf", &x[i]);


        printf("\n Input data 출력 ");

        for (i = 1; i <= n; i = i + 1){

                printf("\n");

                for (j = 1; j <= n; j = j + 1){

                        printf("%10.5f", a[i][j]);

                }

                printf("%10.5f", b[i]);

        }

printf("\n");

        Jacobi();


        printf("\n 계산결과 출력 ");

        printf("\n iter = %d, 평균 residual = %f ", iter, resavg);

        for (i = 1; i <= n; i = i + 1){

                printf("\n x(%d)=%10.6f", i, x[i]);

        }

}


void Jacobi(){

        double ahat[nd][nd], bhat[nd], res[nd];

        double sum, ressum;

        int i, j;

        for (i = 1; i <= n; ++i){

                for (j = 1; j <= n; ++j){

                        ahat[i][j] = a[i][j] / a[i][i];

                }

                bhat[i] = b[i] / a[i][i];

                ahat[i][i] = 0.;

        }


        for (iter = 1; iter <= imax; ++iter){

                for (i = 1; i <= n; i++){

                        x[i] = bhat[i];

                        for (j = 1; j <= n; ++j){

                                x[i] = x[i] - ahat[i][j] * x[j];

                        }

                }

                ressum = 0;

for (i = 1; i <= n; ++i){

                        sum = x[i] - bhat[i];

                        for (j = 1; j <= n; ++j){

                                sum = sum + ahat[i][j] * x[j];

                        }

                        res[i] = fabs(sum);

                        ressum = ressum + res[i];

                }

                resavg = ressum / n;

                if (resavg < epsilon){

                        goto END;

                }

        }

                printf("\n 수렴하지 않는다 ");

        END:;

}



'언어 > Coding' 카테고리의 다른 글

역 거듭제곱  (0) 2015.12.22
Lagrange  (0) 2015.12.22
Jacobi  (0) 2015.12.21
GaussN  (0) 2015.12.21
Least(최소 제곱법)  (0) 2015.12.21
Posted by 知彼知己百戰不殆
,