Metodat numerike

[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();    
    
}