# Shkenca > Informatikë dhe Internet > Arti i programimit >  Metodat numerike

## suli

Pershendetje!!
kerkoj ndihme ne shkrimin e nje programi ne C per Metoden e Jacobit per te gjetur "vlerat e veta te matricave". Ku matrica eshte simetrike dhe me vlera reale.
do ju isha sh mirenjohes
flm - suli

----------


## qoska

Me fal qe po jua them po pse sprovoni me perpara te kerkoni njehere ne google, yahoo per te gjetur keto programe.
Per kete ne fjale ja kam gjetur nje personi para ca kohesh dhe ekzistonin versione me shumice.

----------


## suli

flm per keshillen....

----------


## Ardi_Pg_ID

Gjithashtu te sygjeroj googel po perseri ate script qe do gjesh e googel perdore si baze per te studjuar menyren edhe vijueshmerine e programit te rekomandoj edhe te sygjeroj qe te shkruash kodin tend personal do ndjehesh shume here me mire
Ardi

----------


## suli

nuk eshte me aktuale...
desha vetem te thoja se jam "i ri" ne kete forum, kur shkruajta ketu u nisa me besimin se ndihma do ishte me konkrete se berja moral dhe ca "keshilla"...
anyway

----------


## Ardi_Pg_ID

Ndihma kalon me konkrete ... 
posto kodin qe ke shkruajtur duke ardhur me pyetje me konkrete edhe ne do te te ndihmojme ta permiresosh edhe ta regullosh ate kod 
Ardi

----------


## suli

Ok.

problemi qendron per ta cuar matricen simetrike X->Ak=diag(lambdai);
A1=X;
Ak1 eshte e njejte me matricen Ak, por ndryshon ne elementet
Ak1.App=Ak.App+el.t*Ak.Apq;
Ak1.Aqq=Ak.Aqq-el.t*Ak.Apq;
Ak1.A[j][p]=Ak1.A[p][j]=Ak.A[p][j] - el.s*(Ak.A[q][j] + el.r*Ak.A[p][j]);
Ak1.A[j][q]=Ak1.A[q][j]=Ak.A[q][j] + el.s*(Ak.A[p][j] - el.r*Ak.A[q][j]);
per j!=p
ku t dhe r jane te dhena me formulat: t=c/s r=s/1+c
ku c dhe s jane vlerat e rrotulimeve givensa  me kend fi;
kendi fi llogaritet ne baze elementit me te madh Apq (pervec diagonales) te matrices Ak
ndryshimet nga matrica Ak na Ak1 vazhdojne, deri sa elementet e matrices (pervec diagonales) te jene me te vogla se nje vlere eps.

problemi eshte se ky program nuk punon mire.elementet nuk zerohen.
do te isha mirenjohes po me ndihmove ketu.
te fala suli



```

struct Matrix{
    float A[MAX_SIZE][MAX_SIZE];
    int p;
    int q;
    float Apq;
    float App;
    float Aqq;
};
struct elem{
    float t;
    float r;
    float s;
};


//METODA JACOBIEGO

//funkcje potrzebne

/*
/**
1. funksioni per te gjetru vleren me te madhe te mactices
.kthen kete element, kooordinatat p, dhe q dhe vleren A[p][p] dhe A[q][q]
*/ 

Matrix MaxElement( float A[MAX_SIZE][MAX_SIZE], int n)
{
    Matrix wynik;
    wynik.Apq=0.0;
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++){
            if(A[i][j]>wynik.Apq && i!=j){
                wynik.Apq=A[i][j];
                wynik.p=i;
                wynik.q=j;
                wynik.App=A[i][i];
                wynik.Aqq=A[j][j];
                
            }
            wynik.A[i][j]=A[i][j];
        }
    return wynik;
}

//----------------------------------------
/*
/** funksjoni kthen vlrene t, r dhe s 
*/
elem DifferencElement(Matrix T, int n){
    
    //obliczenie kata
    float fi=0;
    if(T.App-T.Aqq==0){
        if(T.Apq > 0) fi = (float)PI/4;
        else if(T.Apq < 0) fi = (float) -PI/4;
    }else
        fi=0.5*atanf(2*T.Apq/(T.App-T.Aqq));

    //obliczenie wartosci t i r z korymi roznia sie macierz Ak+1
    elem wynik;
    wynik.t=(float)cos(fi)/sin(fi);
    wynik.r=(float)sin(fi)/(1+cos(fi));
    wynik.s=sin(fi);
    return wynik;
}

//----------------------------------------
/*
/**

*/
void overwrite(float A[MAX_SIZE][MAX_SIZE],float B[MAX_SIZE][MAX_SIZE],int n){
    for (int i=0;i<n;i++)
        for(int j=0;j<n;j++)
            A[i][j]=B[i][j];
}

//--------------------------------------
/*
/**
*/
bool IsEnd(Matrix Ak, int n){
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++){
            if(Ak.A[i][j] < eps && i!=j)
                return true;
            else
                return false;
        }
}

//--------------------------------------
/*
/**
funk kthen  Ak+1 ne baze te Ak
*/
Matrix Diagonalizacja(Matrix Ak, int n){
    Matrix Ak1;
    Ak=MaxElement(Ak.A,n);//znajdziemy najwiekszy element w macierzy C (oraz informacje
    int p=Ak.p;
    int q=Ak.q;
    elem el;
    el= DifferencElement(Ak,n);
    //float Ak1[MAX_SIZE][MAX_SIZE];
    overwrite(Ak1.A,Ak.A,n);
    Ak1.App=Ak.App+el.t*Ak.Apq;
    Ak1.Aqq=Ak.Aqq-el.t*Ak.Apq;
    for(int j=0;j<n;j++){
        if(p!=j){
            Ak1.A[j][p]=Ak1.A[p][j]=Ak.A[p][j] - el.s*(Ak.A[q][j] + el.r*Ak.A[p][j]);
            Ak1.A[j][q]=Ak1.A[q][j]=Ak.A[q][j] + el.s*(Ak.A[p][j] - el.r*Ak.A[q][j]);
        }
    }
    return Ak1;
}
//--------------------------------------
/*
/**
Metoda Jacobiego 
*/
void Jacobiego(float C[MAX_SIZE][MAX_SIZE], int n)
{
    Matrix A;
    overwrite(A.A,C,n); //przepisujemy macierz C na A
    elem wart;
    wart=DifferencElement(A,n);
    //sprawdzenie czy jest warunek stopu oraz interacja na macierz Ak+1
    while(IsEnd(A,n)==false){
        A=Diagonalizacja(A,n);
        
        printf("\n\n\n");
    }
    //return A;
} 



```

----------


## Alket123

ky nuk eshte forum matematike. Mburren moderatoret eshte reklame per forumin. Ne mezi sa kemi kaluar klasen.

ke tre shembuj nqs te ndihmojne sadopak.

----------


## Alket123

ose ka nje rruge me te thjeshte. Krijohe ne Matlab, pastaj exportohe/konvertoje ne C

help
http://www.mathworks.com/support/sol...lution=1-18L04

----------


## suli

flm prof.
e vleresoj ndihmen
respekte - suli

----------


## Alket123

Shpjegoje njehere me shume, shkruaje anglisht te kuptoj me mire. Ne Matlab me duket e lehte.

----------


## suli

Pak a shume idea  e metodes se Jacobit eshte qe matrica simetrike C te kthehet ne matrice diagonale. teorikisht kjo zgjat ne infinit, kurse ne ketu do perfundoje kur elementet e matrices ( sigurisht pervec diagonales) te jene me te vogla se nje konstane eps > 0 e deklaruar. E gjithe kjo behet ne menyre interaction.( dmth matrica A(k+1) eshte me afer matrces diagonale se A(k). Si fillim A(1)=C;
Ne matricen A(k) gjejme elementin me te madh (kete element nuk e kerkojme ne diagonale) supozojme se ky eshte elementi A[p][q] ( p, q - indekset). Tani per te formuar matricen A(k+1) perdorim rrotullimet e Givensit me kendin fi.ky kend llogaritet ne baze te formules: tg(2*fi)=2*A[p][q] / (A[p][p] - A[q][q])
nqs A[p][p] - A[q][q] ==0, atehere si kend fi merret +PI/4 ose -PI/4 ne varesi te sgn(A[p][q]) duke njohur kendin fi, gjejme c:=cos(fi) dhe s:=sin(fi) dhe me pas definiojme t:=c/s dhe r:=s/1+c
Matrice A(k+1) ndryshon nga matrica A(k) ne keto elemente:
A(k+1)[p][p] = A(k)[p][p] + t*A(k)[p][q]
A(k+1)[q][q] = A(k)[q][q] -  t*A(k)[p][q]
A(k+1)[j][p] = A(k+1)[p][j] = A(k)[p][j] - s*(A(k)[q][j] + r*A(k)[p][j])   per j=1...n, j!=p
A(k+1)[j][q] = A(k+1)[q][j] = A(k)[q][j] + s*(A(k)[p][j] - r*A(k)[q][j])   per j=1...n, j!=q

ketu eshte dicka qe matrica eshte simetrike, kshu qe kete mund ta perdorim vetem per njeren pjese (psh trekendeshin e siperm)

kjo pak a shume eshte idea

----------


## IlirDeda

Suli, ketu ke kodin per metoden e Jakobit.

http://library.lanl.gov/numerical/bookcpdf/c11-1.pdf 

Provoje njehere nese te ben.

Ai kodi qe ke derguar ketu nuk kompilohet. Me duket se nuk e ke derguar te plote. Nese do qe ta shikojme per gabime do beje mire ta dergoje te gjithin.

----------


## suli

faleminderit Ilir.
e vertete qe kodi nuk eshte i teri, sepse Metoda Jacobi eshte vetem nje pjese e programit ne teresi. puna eshte se versioni i kompiluar eshte sh. i pakuptueshem, mgjth do me duhet ta rregulloj pak rremujen aty, dhe do e vendos ne forum.
flm per linkun
te fala suli

----------


## Alket123

Ku jeni programatore shqiptare ne C ose Matlab? Une ndihmen time e dhashe  me shume sec duhet, me shume sec me takonte, biles po e jap prape duke tundur nga gjumi kete teme per te e vene ne liste te paren.

po shtoj disa tituj te rendesishem ne lidhje me kete teme. Ndofta ne nivel astronomik. Te ju zgjeroj horizontin.

Nonphysical sheet for perturbed Jacobian matrices
http://www.ams.org/mathscinet-getitem?mr=95k:47015

The frequency theorem for the case in which the state space and the control space are Hilbert spaces, and its application in certain problems in the synthesis of optimal control.
http://www.ams.org/mathscinet-getitem?mr=57:13636

Hilbert jane ne Matlab ne menyre fantastike. Do jap nje shebull me vone.

Ne C convertohet nepermjet Matlabit. Ne Matlab nuk ekziston Jacob direkt.
Ndoshta mund te jete e pavlere por "help numjac" ne Matlab

Dikush ishte i interesuar dikur. ne kete teme tregohet pak ndihme ne Matlab per dembelet te cilet nuk zgjidhin ushtrimet vete
si:

x = sym('x')
>> x = sym('x');
>> diff(x^3)
ans =
3*x^2

me shume ne link, tema e fundit.
http://www.forumishqiptar.com/showthread.php?t=46458

----------


## suli

I gjithe programi eshte ky.
nuk jam i sigurte nese eshte i sakte ( per pjesen qe i perket algoritmit te Jacobit, se pj tj eshte ne rregull)
ky si program gjen vlerat e veta ne rastin e ekuacionit
Ax = λBx, ku Matrica A eshte simetrke, dhe B eshte simetrike dhe me vlera pozitive te ndryshme nga zero. (rasti i zakonshem eshte kur B=I)
pak a shume si punon.
pasi jepen vlerat perkatese te matricave, programi matricen B e kthen ne formen 
B=LT*L, ne baze te algoritmit te cholesky - Banachiewicz, ku L eshte matrice e trekendeshit te siperm, dhe LT eshte matrice e transponuar.
me pas per matricen C=L-1*A*(LT)-1 ku (-1) tregon qe eshte matrice e anasjellte, perdor algorytmin e Jacobit.
problemi eshte se per vlera epsilon akoma me te vogla, sikur nuk "zerohen" te gjitha vlerat e matrices.
mgjth do konsultohem akoma, dhe do perpiqem ta permisoj. besoj se se shpejti do kem versionin perfundimtar, qe ndoshta ne te ardhme do ti hyje ne pune ndojerit

pershendetje suli



```

// projekt.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#define MAX_SIZE 20
#define PI 3.14159265358979
#define eps 0.01

//----------------------------------------
struct Matrix{
    float A[MAX_SIZE][MAX_SIZE];
    int p;
    int q;
    float Apq;
    float App;
    float Aqq;
    int epsilon;
};
struct elem{
    float t;
    float r;
    float s;
};
//----------------------------------------

//----------------------------------------
void wyswietlanie_ok(float A[MAX_SIZE][MAX_SIZE], int n)
{
    for(int i=0;i<n;i++){
            for(int j=0;j<n;j++)
                printf("%f\t",A[i][j]);
                printf("\n");
        }
    printf("\n\n");
}
//----------------------------------------
void wyswietlanie(float A[MAX_SIZE][MAX_SIZE], int n)
{
    for(int i=0;i<n;i++){
            for(int j=0;j<n;j++)
                printf("%.2f\t",A[i][j]);
                printf("\n");
        }
    printf("\n\n");
}
//----------------------------------------
Matrix mnozenie(float A[MAX_SIZE][MAX_SIZE],float B[MAX_SIZE][MAX_SIZE], int n){
    Matrix C;
    
    for(int m=0;m<n;m++)
        for(int l=0;l<n;l++)
            C.A[m][l]=0.0;
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++)
            for(int k=0;k<n;k++)
                C.A[i][j]+=A[i][k]*B[k][j];
    return C;
}

//----------------------------------------

void zakoncz(){
    printf("zakoncz\n");
    exit(0);
}



//METODA JACOBIEGO

//funkcje potrzebne

/*
/**
1. funkcja do znalezienia najwiekszego elementu w macierzy;
zwraca ten element, wsporzedne p oraz q, wartosc A[p][p] i A[q][q]
*/ 

Matrix MaxElement( float A[MAX_SIZE][MAX_SIZE], int n)
{
    Matrix wynik;
    wynik.Apq=0.0;
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++){
            if(fabs(A[i][j])>wynik.Apq && i!=j){
                wynik.Apq=A[i][j];
                wynik.p=i;
                wynik.q=j;
                wynik.App=A[i][i];
                wynik.Aqq=A[j][j];                
            }
            wynik.A[i][j]=A[i][j];
        }

        printf("\nky eshte max = %f\n\n",wynik.Apq);
    return wynik;
}

//----------------------------------------
/*
/** funkcja ta zwraca wartosci t i r potrzebne aby wykonywac
nastepne kroki interacji
*/
elem DifferencElement(Matrix T, int n){
    
    //obliczenie kata
    float fi=0;
    if(T.App-T.Aqq==0){
        if(T.Apq > 0) fi = (float)PI/4;
        else if(T.Apq < 0) fi = (float) -PI/4;
    }else
        fi=0.5*atanf(2*T.Apq/(T.App-T.Aqq));

    //obliczenie wartosci t i r z korymi roznia sie macierz Ak+1
    elem wynik;
    wynik.t=(float)cos(fi)/sin(fi);
    wynik.r=(float)sin(fi)/(1+cos(fi));
    wynik.s=sin(fi);
    return wynik;
}

//----------------------------------------
/*
/**
Funkcja przepisujaca macierz B na macierz A
*/

void overwrite(float A[MAX_SIZE][MAX_SIZE],float B[MAX_SIZE][MAX_SIZE],int n){
    for (int i=0;i<n;i++)
        for(int j=0;j<n;j++)
            A[i][j]=B[i][j];
}


//--------------------------------------
/*
/**
funkcja ktora zwraca macierz Ak+1 od macierzy Ak
*/
Matrix Diagonalizacja(Matrix Ak, int n){
    Matrix Ak1;
    Ak=MaxElement(Ak.A,n);//znajdziemy najwiekszy element w macierzy C (oraz informacje
    int p=Ak.p;
    int q=Ak.q;
    elem el;
    el= DifferencElement(Ak,n);
    printf("p=%i, q=%i, Apq=%f, App%f, Aqq=%f\n\n",Ak.p,Ak.q,Ak.Apq,Ak.App,Ak.Aqq);
    printf("r=%f, t=%f, s=%f\n\n",el.r,el.t,el.t);
    //float Ak1[MAX_SIZE][MAX_SIZE];
    overwrite(Ak1.A,Ak.A,n);
    Ak1.App=Ak.App+el.t*Ak.Apq;
    Ak1.Aqq=Ak.Aqq-el.t*Ak.Apq;
    for(int j=0;j<n;j++){
        if(p!=j){
            Ak1.A[j][p]=Ak1.A[p][j]=Ak.A[p][j] - el.s*(Ak.A[q][j] + el.r*Ak.A[p][j]);
            Ak1.A[j][q]=Ak1.A[q][j]=Ak.A[q][j] + el.s*(Ak.A[p][j] - el.r*Ak.A[q][j]);
        }
    }
    wyswietlanie_ok(Ak1.A,n);
    Ak1.epsilon=0;
    
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++)
            if(Ak.A[i][j]<eps && i!=j){
                Ak1.epsilon++;
            }
            printf("\n\neps=%i\n",Ak1.epsilon);
    return Ak1;
}
//--------------------------------------
/*
/**
Metoda Jacobiego 
*/
Matrix Jacobiego(float C[MAX_SIZE][MAX_SIZE], int n)
{
    Matrix A,Ak;
    overwrite(A.A,C,n); //przepisujemy macierz C na A
    Ak=Diagonalizacja(A,n);
    while(1)
    {
        if(Ak.epsilon<n*n-n)
            Ak=Diagonalizacja(Ak,n);
        else{
            wyswietlanie_ok(Ak.A,n);
            return Ak;
        }
    }
}
//--------------------------------------
/*
/**
fukncja ktora zwraca macierz odwrotna macierzy trojkatnej dolnej L
*/

Matrix Odwrotna( float A[MAX_SIZE][MAX_SIZE], int n){
    Matrix X;
    for(int m=0;m<n;m++)
        for(int l=0;l<n;l++)
            X.A[m][l]=0.0;
    float Z[MAX_SIZE][MAX_SIZE];
    for(int i=0;i<n;i++){
        X.A[i][i]=(float)1/A[i][i];
        for(int k=i+1;k<n;k++)
            for(int j=i;j<k;j++)
                X.A[k][i]=(float)-A[k][j]*X.A[j][i]/A[k][k];
    }
    return X;
}


//----------------------------------------
void obliczenia(){
    float A[MAX_SIZE][MAX_SIZE]={0},B[MAX_SIZE][MAX_SIZE];
int n,i,j;
    printf("obliczenia\n\n");
    printf("podaj rzad macierzy ");
    scanf("%i",&n);
    if(n<=0){
        printf("blad, rzad musi byc dodatny\n");
        exit(0);
    }

    printf("\ndane dla macierzy symetrycznej A\n\n");
    for (i=0;i<n;i++)
        for (j=i;j<n;j++)
        {
            printf("A[%d][%d]: ",i+1,j+1);
            scanf("%f",&A[i][j]);
        }
        printf("\n");
        printf("\n");
    for(i=n-1;i>=0;i--)
        for(j=n-1;j>=i;j--)
            A[j][i]=A[i][j];
    
    printf("\ndane dla macierzy symetrycznej, dodatnio okreslonej B\n\n");
    for (i=0;i<n;i++)
        for (j=i;j<n;j++)
        {
            printf("B[%d][%d]: ",i+1,j+1);
            scanf("%f",&B[i][j]);
            if(B[i][j]<=0)
            {
                printf("Macierz B jest musi byc dodatnio okreslona\n");
                printf("B[%d][%d]: ",i+1,j+1);
                scanf("%f",&B[i][j]);
                if(B[i][j]<=0)
                    exit(0);
            }
        }
    for(i=n-1;i>=0;i--)
        for(j=n-1;j>=i;j--)
            B[j][i]=B[i][j];

        printf("\n\nMacierz A\n");
        wyswietlanie(A,n);
        printf("\n\nMacierz B\n");
        wyswietlanie(B,n);
        printf("\n\n");

        //algorytm choleskyego - Banachiewicza
        float L[MAX_SIZE][MAX_SIZE]={0.0};
        for(int k=0;k<n;k++)
        {
            if(B[k][k]<0.0){
                printf("%f < 0.0\n",L[k][k]);
            //    exit(0);
            }
            else
            {
                L[k][k]=sqrt((float)B[k][k]);
                for(int s=k+1;s<n;s++)
                    if(L[k][k]==0.0)printf("dzielenie przez 0\n");
                    else
                    L[s][k]=(float)B[s][k]/L[k][k];
                for(j=k+1;j<n;j++)
                    for(i=j;i<n;i++)
                        B[i][j] = fabs(B[i][j] - L[i][k]*L[j][k]); //ketu kam marre fabs-duhet?? pyete
            }
        }
        
        L[n][n]=sqrt((float)B[n][n]);
        //----------------------------
        printf("Macierz trojdiagonalna L\n");
        wyswietlanie(L,n);
        float LT[MAX_SIZE][MAX_SIZE];
        for(i=0;i<n;i++){
            for(j=0;j<n;j++)
                LT[i][j]=L[j][i];
        }
        printf("\n\nMAcierz LT\n");
        wyswietlanie(LT,n);
        Matrix S,SS,I,II; //S=(L)-1; SS=(LT)-1
        //TW. (AT)-1=(A-1)T
        S=Odwrotna(L,n);
        printf("\nodwrotna matricy L\n");
        wyswietlanie_ok(S.A,n);
        
        for(i=0;i<n;i++){
            for(j=0;j<n;j++)
                SS.A[i][j]=S.A[j][i];
        }

        wyswietlanie_ok(SS.A,n);

        
        tmp=mnozenie(S.A,A,n);
        wyswietlanie(tmp.A,n);
        C=mnozenie(tmp.A,SS.A,n);
        printf("macierz C dla ktorej bedzie\n");
        printf("zastosowana metoda Jacobiego\n\n");
        wyswietlanie_ok(C.A,n);
        

        temp=Jacobiego(C.A,n);
        printf("\n\nWartosci walsne\n");
        for(i=0;i<n;i++)
            printf("%f\n",temp.A[i][i]);

        printf("\n\n\n");
        
    menu();    
    
}

void main(int argc, _TCHAR* argv[])
{
    
     obliczenia();    
    
} 



```

----------

