ot_systeme.cpp

Aller à la documentation de ce fichier.

//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//  MAGiC
//####//  Jean Christophe Cuilliere et Vincent FRANCOIS
//####//  Departement de Genie Mecanique - UQTR
//####//------------------------------------------------------------
//####//  MAGIC est un projet de recherche de l equipe ERICCA
//####//  du departement de genie mecanique de l Universite du Quebec a Trois Rivieres
//####//  http://www.uqtr.ca/ericca
//####//  http://www.uqtr.ca/
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//
//####//       ot_systeme.cpp
//####//
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//    COPYRIGHT 2000-2024
//####//  jeu 13 jun 2024 11:53:59 EDT
//####//------------------------------------------------------------
//####//------------------------------------------------------------
 
 
 
#pragma hdrstop
#include <stdlib.h>
#include <math.h>
#include "ot_systeme.h"
#include "ot_mathematique.h"
 
 
OT_SYSTEME::OT_SYSTEME()
{
}
 
 
 
OT_SYSTEME::~OT_SYSTEME()
{
}
 
 
int OT_SYSTEME::resoud_gauss(double *a,double* b,double *x,int n)
{
    double *aaug=new double[(n+1)*n];
#define a(i,j) (*(a+(i)*(n)+j))
#define aaug(i,j) (*(aaug+(i)*(n+1)+j))
    for (int i=0;i<n;i++)
        for (int j=0;j<n;j++)
            aaug(i,j)=a(i,j);
    for (int k=0;k<n;k++)
        aaug(k,n)=b[k];
    for (int j=0;j<n-1;j++)
    {
        // recherche le pivot
        double pvt=fabs(aaug(j,j));
        int ligne_pivot_debut=j;
        int ligne_pivot=j;
 
        for (int i=j+1;i<n;i++)
            if (fabs(aaug(i,j))>pvt)
            {
                pvt=fabs(aaug(i,j));
                ligne_pivot=i;
            }
        if (ligne_pivot!=ligne_pivot_debut) // trouver meilleur
        {
            for (int col=j;col<n+1;col++)
            {
                double temp=aaug(ligne_pivot_debut,col);
                aaug(ligne_pivot_debut,col)=aaug(ligne_pivot,col);
                aaug(ligne_pivot,col)=temp;
            }
        }
        // j ai mon pivot max
        if (fabs(aaug(j,j))< 0.00000000000001) {
            delete [] aaug;
            return 0;
        };
        // j ai mon pivot max non nul
        for (int ii=j+1;ii<n;ii++) aaug(ii,j)=aaug(ii,j)/aaug(j,j);
 
        for (int iii=j+1;iii<n;iii++)
            for (int k=j+1;k<n+1;k++)
                aaug(iii,k)=aaug(iii,k)-aaug(iii,j)*aaug(j,k);
    }
 
    if (fabs(aaug(n-1,n-1))<0.00000000000001) {
        delete [] aaug;
        return 0;
    };
    x[n-1]=aaug(n-1,n)/aaug(n-1,n-1);
    for (int j=n-2;j>-1;j--)
    {
        x[j]=aaug(j,n);
        for (int k=n-1;k>j;k--)
            x[j]=x[j]-aaug(j,k)*x[k];
        x[j]=x[j]/aaug(j,j);
    }
 
 
    delete [] aaug;
    return 1;
#undef a // enlever le  (i,j) pour enlever le warning
#undef aaug
}
 
int OT_SYSTEME::resoud_QR(double *m, double *b, double *x, int dim1,int dim2)
{
    double *mat=new double[dim1*(dim2+1)];
    double *r=new double[dim2];
#define mat(i,j) (*(mat+(i)*(dim2+1)+(j)))
#define m(i,j) (*(m+(i)*dim2+(j)))
    for (int i=0;i<dim1;i++)
    {
        for (int j=0;j<dim2;j++)
            mat(i,j)=m(i,j);
        mat(i,dim2)=b[i];
    }
    for (int k=0;k<dim2;k++)
    {
        double somme=0.;
        for (int i=k;i<dim1;i++)
            somme=somme+mat(i,k)*mat(i,k);
        double alpha=sqrt(somme);
        double pivot=mat(k,k);
        if (pivot>0) alpha=-alpha;
        r[k]=alpha;
        if (OPERATEUR::egal(alpha,0.,0.000001))
        {
            delete [] mat;
            delete [] r;
            return 0;
        }
        double beta=somme-pivot*alpha;
        pivot=pivot-alpha;
        mat(k,k)=pivot;
        for (int j=k+1;j<dim2+1;j++)
        {
            double somme=0.;
            for (int i=k;i<dim1;i++)
                somme=somme+mat(i,k)*mat(i,j);
            double gamma=somme/beta;
            for (int i=k;i<dim1;i++)
                mat(i,j)=mat(i,j)-gamma*mat(i,k);
        }
    }
    x[dim2-1]=mat(dim2-1,dim2)/r[dim2-1];
    for (int j=dim2-2;j>-1;j--)
    {
        x[j]=mat(j,dim2);
        for (int k=dim2-1;k>j;k--)
            x[j]=x[j]-mat(j,k)*x[k];
        x[j]=x[j]/r[j];
    }
 
 
 
#undef mat
    delete [] mat;
    delete [] r;
    return 1;
}
 
 
 
 
 
 
int OT_SYSTEME::decompose_LU(double *a,double* res,int dim)
{
#define a(i,j) (*(a+(i)*(dim)+(j)))
#define res(i,j) (*(res+(i)*(dim)+(j)))
 
    if (OPERATEUR::egal(a(0,0),0.,0.00001)==1) return 0;
    for (int i=0;i<dim;i++)
        res(i,0)=a(i,0);
    for (int j=1;j<dim;j++)
        res(0,j)=a(0,j)/res(0,0);
 
 
    for (int j=1;j<dim;j++)
    {
        for (int i=j;i<dim;i++)
        {
            double sum=0.;
            for (int k=0;k<j;k++)
                sum=sum+res(i,k)*res(k,j);
            res(i,j)=a(i,j)-sum;
        }
        for (int i=j+1;i<dim;i++)
        {
            double sum=0.;
            if (OPERATEUR::egal(res(j,j),0.,0.00001)==1) return 0;
            for (int k=0;k<j;k++)
                sum=sum+res(j,k)*res(k,i);
            res(j,i)=(a(j,i)-sum)/res(j,j);
        }
 
    }
    return 1;
 
#undef a
#undef res
}
 
 
int OT_SYSTEME::resoud_LU(double *matlu,double *b,double *x,int dim)
{
#define matlu(i,j) (*(matlu+(i)*(dim)+(j)))
 
    double* bb=new double[dim];
    bb[0]=b[0]/matlu(0,0);
    for (int i=1;i<dim;i++)
    {
        double sum=0.0;
        for (int k=0;k<i;k++)
            sum=sum+matlu(i,k)*bb[k];
        bb[i]=(b[i]-sum)/matlu(i,i);
    }
 
    x[dim-1]=bb[dim-1];
    for (int j=dim-2;j>=0;j--)
    {
        double sum=0.0;
        for (int k=j+1;k<dim;k++)
            sum=sum+matlu(j,k)*x[k];
        x[j]=bb[j]-sum;
 
    }
 
#undef lu
 
    delete [] bb;
    return 1;
}
 
 
 
 
 
 
void OT_SYSTEME::get_inv(double *mat1,double *mat2,int n,int *ierr)
{
#define mat(i,j) (*(mataug+(j)*n+i))
#define mat1(i,j) (*(mat1+(j)*n+i))
#define mat2(i,j) (*(mat2+(j)*n+i))
 
    double *mataug=(double*)calloc(2*n*n,sizeof(double));
 
    for (int i=0;i<n;i++)
    {
        for (int j=0;j<n;j++)
            mat(i,j)=mat1(i,j);
        mat(i,i+n)=1.;
    }
 
 
    *ierr=0;
 
    for (int j=0;j<n;j++)
    {
        double pivot=fabs(mat(j,j));
        int ligne_pivot=j;
        for (int i=j+1;i<n;i++)
        {
            if (fabs(mat(i,j))>pivot)
            {
                pivot=fabs(mat(i,j));
                ligne_pivot=i;
            }
        }
        if (pivot<0.000000001) {
            *ierr=1;
            free(mataug);
            return;
        }
        if (ligne_pivot!=j)
        {
            for (int i=j;i<2*n;i++)
            {
                double temp=mat(j,i);
                mat(j,i)=mat(ligne_pivot,i);
                mat(ligne_pivot,i)=temp;
            }
        }
        for (int i=j+1;i<n;i++)
            mat(i,j)=-mat(i,j)/mat(j,j);
        for (int i=j+1;i<n;i++)
            for (int k=j+1;k<2*n;k++)
                mat(i,k)=mat(i,k)+mat(i,j)*mat(j,k);
    }
 
    for (int i=0;i<n;i++)
    {
        mat2(n-1,i)=mat(n-1,n+i)/mat(n-1,n-1);
        for (int j=n-2;j>-1;j--)
        {
            mat2(j,i)=mat(j,n+i);
            for (int k=n-1;k>j;k--)
                mat2(j,i)=mat2(j,i)-mat2(k,i)*mat(j,k);
            mat2(j,i)=mat2(j,i)/mat(j,j);
        }
    }
 
    free(mataug);
#undef mat
#undef mat1
#undef mat2
}