Don't like ads? PRO users don't see any ads ;-)
Public Pastes
Guest

Untitled

By: a guest on Jun 1st, 2012  |  syntax: None  |  size: 7.92 KB  |  hits: 12  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1. #include <iostream>
  2. #include <math.h>
  3. #include <iomanip>
  4. #include<ctime>
  5.  
  6.  
  7.  
  8. using namespace std;
  9. void Cholesky(double ** A, int taille,double * B,double * sk);
  10. void Affichage_matrice(double ** A, int taille);
  11. int main()
  12. {    clock_t Time;
  13.  
  14.     int nbr_coef=3;
  15.     double vk[3]={0,0,0}; // Vk 3 coef dans le cercle, 4 pour l'amortisseur
  16.     double sk[3]={0,0,0};
  17.      double* W=new double[10];
  18.     double W1[3];
  19.      int nbr_points=6; // 6 dans le cas du cercle, 10 sinon!
  20.  
  21.  
  22.  
  23.  
  24.  
  25.  
  26.  
  27.  
  28.     ////********************************* Declaration dynamique des matrices *********************************////
  29.  
  30.     double ** A;
  31.     A = new double* [nbr_coef];
  32.     for ( int i=0 ; i < nbr_coef ; i++)
  33.     A[i] = new double [nbr_coef];
  34.  
  35.  
  36.  
  37.  
  38.    for(int i=0;i<nbr_coef;i++ )
  39.     {
  40.         W1[i]=0;
  41.  
  42.         for(int j=0;j<nbr_coef;j++)
  43.             A[i][j]=0;
  44.  
  45.  
  46.     }
  47.  
  48.     double B[6][3]={NULL};
  49.     double BT[3][6]={NULL};
  50.  
  51.  
  52.    //Cholesky(A,nbr_coef,W1,sk);
  53.  
  54.  
  55.      double X[6]={1,2,3,5,7,9};
  56.      double Y[6]={7,6,7,8,7,5};
  57.  
  58.  
  59.  
  60.       for(int j=0;j<nbr_points;j++)
  61.      {
  62.  
  63.           vk[1]+=(1.0/nbr_points)*X[j];
  64.           vk[2]+=(1.0/nbr_points)*Y[j];
  65.           vk[3]+=X[j]*X[j]+Y[j]*Y[j];
  66.       }
  67.        vk[3]=sqrt(vk[3])*(1.0/nbr_points);
  68.  
  69.  
  70.     //double temps[10]={0.1,0.3,0.7,1.2,1.6,2.2,2.7,3.1,3.5,3.9};
  71.     //double yi[10]={0.558,0.569,0.176,-0.207,-0.133,0.132,0.055,-0.09,-0.069,0.027};
  72.     double norme=2;
  73.  
  74. int z=0;
  75.     while(z<30)
  76.     {z++;
  77.         norme=0;
  78.         for(int i=0;i<nbr_coef;i++ )
  79.         {
  80.             W1[i]=0;
  81.  
  82.             for(int j=0;j<nbr_coef;j++)
  83.                 A[i][j]=0;
  84.  
  85.  
  86.         }
  87.  
  88.         ////********************************* Calcul de la jacobienne pour l'amortisseur *********************************////
  89. /*
  90.     for(int i=0;i<nbr_coef;i++)
  91.     {
  92.  
  93.  
  94.         B[i][0]=exp(-vk[1]*temps[i])*sin(vk[2]*temps[i]+vk[3]);
  95.         B[i][1]=-vk[0]*temps[i]*exp(-vk[1]*temps[i])*sin(vk[2]*temps[i]+vk[3]);
  96.         B[i][2]=temps[i]*vk[0]*exp(-vk[1]*temps[i])*cos(vk[2]*temps[i]+vk[3]);
  97.         B[i][3]=vk[0]*exp(-(vk[1]*temps[i]))*cos(vk[2]*temps[i]+vk[3]);
  98.         }
  99.  
  100.         ////********************************* Calcul de la Jacobienne pour le cercle *********************************////
  101.  
  102.         for(int i=0;i<nbr_points;i++)
  103.         {
  104.  
  105.  
  106.             B[i][0]=(vk[1]-X[i])/(sqrt((vk[1]-X[i])*(vk[1]-X[i])+(vk[2]-Y[i])*(vk[2]-Y[i])));
  107.             B[i][1]=(vk[2]-X[i])/(sqrt((vk[1]-X[i])*(vk[1]-X[i])+(vk[2]-Y[i])*(vk[2]-Y[i])));
  108.             B[i][2]=-1;
  109.  
  110.             }
  111.  
  112.         for(int i=0;i<nbr_points;i++)
  113.         {
  114.             cout<<endl;
  115.             for(int j=0;j<nbr_coef;j++)
  116.                 cout<<B[i][j]<<"  ";
  117.  
  118.         }
  119.  
  120.  
  121.     ////********************************* Calcul de BT *********************************////
  122.  
  123.     for(int i=0;i<nbr_points;i++)
  124.     {
  125.         for(int j=0;j<nbr_coef;j++)
  126.         {
  127.             BT[j][i]=B[i][j];
  128.  
  129.  
  130.         }
  131.  
  132.     }
  133. ////********************************* Declaration de W amortisseur *********************************////
  134.    /* for(int i=0;i<nbr_coef;i++)
  135.     {
  136.         W[i]=yi[i]-vk[0]*exp(-vk[1]*temps[i])*sin(vk[2]*temps[i]+vk[3]);
  137.  
  138.     }
  139.  
  140.     ////********************************* Declaration de W cercle *********************************////
  141.  
  142.     for(int i=0;i<nbr_points;i++)
  143.         {
  144.  
  145.         W[i]=(sqrt((vk[1]-X[i])*(vk[1]-X[i])+(vk[2]-Y[i])*(vk[2]-Y[i])))-vk[3];
  146.         cout<<" bite:"<<W[i]<<" ";
  147.         }
  148.  
  149.  
  150.  
  151.  
  152.     ////********************************* Calcul de BT W*********************************////
  153.  
  154.     for(int i=0;i<nbr_coef;i++)
  155.     {
  156.  
  157.         for(int j=0;j<nbr_points;j++)
  158.         {
  159.             W1[i]+=BT[i][j]*W[j];
  160.  
  161.         }
  162.     }
  163.  
  164.  
  165.  
  166.  
  167. for(int i=0;i<nbr_coef;i++)
  168. {
  169.     for(int j=0;j<nbr_coef;j++)
  170.     {
  171.         for(int k=0;k<nbr_points;k++)
  172.         {
  173.  
  174.                 A[i][j]+=BT[i][k]*B[k][j];
  175.  
  176.  
  177.         }
  178.  
  179.     }
  180.  
  181.  
  182.  
  183. }
  184.  
  185.  
  186.  
  187.          Cholesky(A,nbr_coef,W1,sk);
  188.  
  189.          for(int i=0;i<nbr_coef;i++)
  190.             {
  191.  
  192.              vk[i]+=sk[i];
  193.  
  194.          cout<<sk[i];
  195.          }
  196.  
  197.  
  198.          ////********************************* Calcul de la norme pour l'amortisseur*********************************////
  199.  
  200.  
  201.          for(int i=0;i<nbr_coef;i++)
  202.          {
  203.              norme+=W1[i]*W1[i];
  204.          }
  205.          norme=sqrt(norme);
  206.  
  207.  
  208.  
  209.  
  210.  
  211.     }
  212.     Time=clock();
  213.     cout<<"norme: "<<norme<<" nombre d'iterations : "<<z<<" en "<<(double) Time*1000.0/CLOCKS_PER_SEC <<" ms";
  214.     for(int i=0;i<nbr_coef;i++)
  215.     {
  216.         cout<<endl<<vk[i];
  217.     }
  218.  
  219.  
  220.     ////********************************* Libération de la mémoire *********************************////
  221.  
  222.  
  223.     for ( int i=0 ; i < nbr_coef ; i++)
  224.     delete A [i];
  225.     delete A;
  226.  
  227.    return 0;
  228. }
  229.  
  230.  ////********************************* Fonction d'affichage *********************************////
  231. void Affichage_matrice(double **A, int taille)
  232. {
  233.     int i,j;
  234.  
  235.  
  236.  
  237.     for(i=0;i<taille;i++)
  238.     {
  239.         for(j=0;j<taille;j++)
  240.         {
  241.             cout<<A[i][j]<<" ";
  242.         }
  243.                   cout<<endl;
  244.     }
  245. cout<<endl;
  246.  
  247. }
  248.  
  249. void Cholesky(double ** A,int taille, double *B,double * sk)
  250. {
  251.  
  252.  
  253.     int n,i,j,k;
  254.  
  255.  
  256.     n=0;
  257.  
  258.     ////********************************* Declaration de A X et B *********************************////
  259.         double L[taille][taille],LT[taille][taille],y[taille],X[taille];;
  260.         for(i=0;i<taille;i++)
  261.         {
  262.             for(j=0;j<taille;j++)
  263.             {
  264.                 L[i][j]=0;
  265.                 LT[i][j]=0;
  266.  
  267.             }
  268.  
  269.         }
  270.  
  271.  
  272.  
  273.  
  274.  
  275.  
  276.     ////********************************* Decomposition de A en LLT *********************************////
  277.         for(i=0;i<taille;i++)
  278.         {
  279.             for(j=0;j<=i;j++)
  280.             {
  281.                 if(i==j)
  282.                 {
  283.                     L[i][j]=A[i][j];
  284.                     if(j>=1)
  285.                     {
  286.                     for(k=0;k<=(j-1);k++)
  287.                     {
  288.                         L[i][j]=L[i][j]-(L[i][k]*L[i][k]);
  289.                     }
  290.                     }
  291.  
  292.                     L[i][j]=sqrt(L[i][j]);
  293.  
  294.  
  295.                 }
  296.                 else
  297.                 {
  298.                     L[i][j]=A[j][i];
  299.  
  300.                     if(j>=1)
  301.                          {
  302.                     for(k=0;k<=(j-1);k++)
  303.                     {
  304.                         L[i][j]=L[i][j]-(L[j][k]*L[i][k]);
  305.  
  306.                     }
  307.                     }
  308.                     L[i][j]/=L[j][j];
  309.  
  310.                 }
  311.             }
  312.         }
  313.  
  314.  
  315.         for(i=0;i<taille;i++)
  316.         {
  317.             for(j=0;j<taille;j++)
  318.             {
  319.                 LT[i][j]=L[j][i];
  320.  
  321.             }
  322.  
  323.         }
  324.     ////********************************* Résolution du système Ly=b *********************************////
  325.         y[0]=B[0]/L[0][0];
  326.  
  327.  
  328.         for(i=1;i<taille;i++)
  329.         {
  330.             y[i]=B[i];
  331.             for(j=0;j<i;j++)
  332.             {
  333.                 y[i]-=(y[j]*L[i][j]);
  334.  
  335.             }
  336.             y[i]/=L[i][i];
  337.  
  338.         }
  339.  
  340.     ////********************************* Résolution du système y=LTX *********************************////
  341.  
  342.         X[taille-1]=y[taille-1]/LT[taille-1][taille-1];
  343.  
  344.         for(i=taille-2;i>=0;i--)
  345.         {
  346.             X[i]=y[i];
  347.             for(j=taille-1;j>i;j--)
  348.             {
  349.                 X[i]+=-X[j]*LT[i][j];
  350.  
  351.             }
  352.             X[i]/=LT[i][i];
  353.  
  354.         }
  355.  
  356.  
  357.     ////********************************* Gestion affichage *********************************////
  358.  
  359.  
  360.        /* for(i=0;i<taille;i++)
  361.         {
  362.             cout<<endl;
  363.             for(j=0;j<taille;j++)
  364.             {
  365.                 cout<<L[i][j]<<" ";
  366.             }
  367.         }
  368.  
  369.             cout<<endl;
  370.             cout<<"X";
  371.             cout<<endl;
  372.             for(i=0;i<taille;i++)
  373.             {
  374.                 cout<<B[i]<<"\t"<<X[i]<<"\t"<<y[i] <<endl;
  375.  
  376.  
  377.             }*/
  378.         for(i=0;i<4;i++)
  379.         {
  380.             sk[i]=X[i];
  381.         }
  382.  
  383.  
  384.  
  385.  
  386. }
create a new version of this paste RAW Paste Data