Discussion :

[the_reward]

Bonjour,

Dans mon projet j'essaye de résoudre un système linéaire de type AX=CY où C est une matrice pour le solveur j'ai choisi le gradient conjugué.

Comme ma matrice elle est pas SDP j'ai multiplié par sa transposée donc je dois résoudre tAAX=tACY.

le problème c que la solution converge dès la 2ème itérations (c louche ) et aussi le code ne marche que pour un temps final=0.01!!!!

svp si quelqu'un pourra m'aider je serai vraiment reconnaissante.

Merci d'avance et bonne fin de journée
Voilà mon code:
le main

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
 
program main
 
  use parametremod
  use matricemod
  use pmvmod
  use gcmod
 
  implicit none
 
  integer                                 :: n, i,k, it,nt
  real*8                                  :: xdeb, xfin, l, abcs, temps, temps_fin, pi,z,w,s,q,sigma,r
  real*8, dimension(:) , allocatable      ::x1,y,res,res1,res2,res3,res0, X, U, v, smb, Un,ll
  real*8, dimension( : , : ), allocatable :: M, M2, tM
 
  !lecture des données
  open(UNIT = 10, FILE = 'parametre')
  read(10,*) xdeb, xfin, n, nt, h, alpha, g, temps_fin
  close(10)
 
  !variables
  l = xfin-xdeb
  dx = l*1.D0/n
  dt = temps_fin/(nt-1)
  temps = 1.D-2
 
  a1 = (alpha*dt)/12.d0
  a2 = (1.d0-(1.d0/alpha))*g
  a = a1*a2
  b1 = (g*dt*dx**2)/(4.d0*h**2)
  b2 = (alpha*dt)/6.d0
  b = b1+(b2*a2)
  c = (alpha*dx)/3.d0
  d1 = (dx**3)/(h**2)
  d2 = (2*alpha*dx)/3.d0
  d = d1+d2
  e = (dt*h)/(4.d0*dx)
  !e = (dt*dx**2)/(4.d0*h)
  !f = (dx**3)/(h**2)
  r = 1.d0/((alpha*h**2*(a2/g))/3.d0)
  !print*,'a1=',a1
  !print*,'a2=',a2
  print*,'a=',a
  !print*,'b1=',b1
  !print*,'b2=',b2
  print*,'b=',b
  print*,'c=',c
  !print*,'d1=',d1
  !print*,'d2=',d2
  print*,'d=',d
  print*,'e=',e
  !print*,'f=',f
  print*,dx
  print*,dt
  print*,l
  print*,n
 
 
  !allocation memoire
 
  allocate(X(2*n), U(2*n), v(2*n), Un(2*n), smb(2*n),x1(2*n),y(2*n),res(2*n),res1(2*n),res2(2*n)&
&,res3(2*n),res0(2*n))
  allocate(M(2*n,2*n), M2(2*n,2*n), tM(2*n,2*n))
 
  !condition initiale
  open(UNIT = 10, FILE = 'solution_ini.dat')
!==========================================================================================
  pi = 4.D0*datan(1.D0)
  !print*,pi
  sigma = 0.15d0
  do i = 1, n
     abcs = (i-1) * dx - l/2.D0
     ! zeta gaussienne
     U(i) = 1.D0 / sqrt(2.D0*pi) / sigma * exp(-.5d0 * abcs**2 / sigma**2)
     ! u gaussienne
     U(n+i) = 0.d0!1.D0 / sqrt(2.D0*pi) / sigma * exp(-.5d0 * abcs**2 / sigma**2)
     write(10,*) abcs, U(i)
  end do
  !==================================================================================
  !do i = 1,n
   !  U(i) = 2*cosh((1.d0/sqrt(r))*(i-1)*dx)
   !  U(n+i) = 0.d0
   !  write(10,*) (i-1)*dx,U(i)
  !enddo
  close(10)
 
  !test si (tAAX,Y) = (AX,AY) pour avoir une matrice SDP
  x1(1:5)=3.d0
  x1(6:n-5)=1.d0
  x1(n-4:2*n-3)=15.d0
  x1(2*n-2:2*n)=-7.d0
  !print*,x1
  !print*,''
 
  call matrice1(M,n)
  !print*,M
 
  res=matmul(M,x1)
  !calcul de (A,X)
  !call pmv1(x1,M,res)
  !print*,res
 
  call transposemat(n,M,tM)
  !print*,tM
 
  !calcul de (tA,RES)
  res1=matmul(tM,res)
  !call pmv2(tM,res,res1)
  y(1:4)=2.d0
  y(5:n-7)=-3.d0
  y(n-6:2*n-8)=5.d0
  y(2*n-7:2*n)=-6.d0
 
  !calcul de (RES1,Y)
  z = dot_product(res1,y)
 
  !calcul de (A,Y)
  !call pmv1(y,M,res2)
  res2=matmul(M,y)
  !calcul de (RES,RES2)
  w = dot_product(res,res2)
  print*,'z=',z
  print*,'w=',w
 
  !si la différence z-w est de l'ordre 10^-14 ma matrice tAA est bien SDP
  s=z-w
  print*,s
 
 
  !résolution du système
  X=0.D0
  Un = U
  do k = 1,nt
     call matrice2(M2,n)
     !print*,M2
     v=matmul(M2,Un)
     !call sm(M2,Un,v)
     !print*,v
     !call pmv2(tM,v,smb)
     smb=matmul(tM,v)
     !print*,smb
     X = gc(M,tM,smb)
     !print*,X
     !call pmv3(M,tM,X,res3)
     res0=matmul(M,X)
     res3=matmul(tM,res0)
     !print*,res3
     !q=tAAX-B
     q=maxval(abs(res3-smb))
     print*,q,maxval(abs(X)),maxval(abs(Un))
     Un=X
     !print*,X
  end do
  !print*,X
  !écriture dans un fichier de la solution finale 
  open(UNIT=12, FILE='solution_fineta.dat')
  do i=1,n
     write(12,*) (i-1)*dx, X(i)
  end do
  close(12)
 
  open(UNIT=12, FILE='solution_finU.dat')
  do i=n+1,2*n
     write(12,*) (i-1)*dx, X(i)
  end do
  close(12)
  deallocate(M,tM,M2,X,U,v,smb,Un)
end program

le grand conj

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
 
module gcmod
contains
 
  function gc(M,tM,s)
    use parametremod
    use matricemod
    !use pmvmod
    implicit none
 
    real*8, dimension(:), intent(in)   :: s
    real*8, dimension(:,:),intent(in)  :: M, tM
    integer                            :: k,nt,n
    real*8, dimension(size(s))         :: gc, x, r,p, q, ancienr, res1,res2,res0,z,w
    real*8                             :: alpha1, beta1,teta,y, err,err1, eps,norme,aa,bb
 
    !initialisation
    x = 0.D0
    q  = 0.D0
    alpha1 = 0.D0
    beta1 = 0.D0
    eps = 1.D-7
    err = 1.D0
    k = 0
    !itérations
    !print*,tM
    !print*,M
    res2=matmul(M,X)
    res1=matmul(tM,res2)
    !call pmv3(M,tM,x,res1)
    !print*,res1
    q = s - res1
    r = q
    do while (err > eps) 
       !call pmv3(M,tM,q,res1)
       res0=matmul(M,q)
       res1=matmul(tM,res0)
       !print*,res1
       alpha1 = dot_product(r,r)/dot_product(res1,q)
       x = x+alpha1*q
       ancienr = r
       r = r -alpha1*res1
       beta1 = dot_product(r,r)/dot_product(ancienr,ancienr)
       q = r +beta1*q
       err1 = dot_product(r,r)
       err = sqrt(err1)
       !print*,err1,err 
       !print*,err
       k = k+1
    end do
    gc = x
    print*,'k=',k
 
  end function gc
 
end module gcmod

Les matrices

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
 
module matricemod
  contains
 
    subroutine matrice1(mat,n)
      use parametremod
      integer, intent(in) :: n
      real*8, dimension(2*n,2*n), intent(out) :: mat
      integer :: i
 
      ! Initialisation
      mat(:,:) = 0
      !print*,size(a,1)
 
      ! bloc 1
      do i = 1,n!size(a,1)/2
         mat(i,i) = 1
      end do
 
      !bloc 2
      mat(n+2,1)=e                           !mat((size(a,1)/2)+2,1)=7
      mat(2*n,1)=-e                        !a(size(a,1),1)=-7
      do i = 2,n-1!(size(a,1)/2)-1
         mat(i+n+1,i) = e                    !a(i+(size(a,1)/2)+1,i) = 7
         mat(i+n-1,i) = -e                     !a((size(a,1)/2)+i-1,i) = -7
      end do
      mat(n+1,n)=e                         !a((size(a,1)/2)+1,size(a,1)/2)=7
      mat(2*n-1,n)=-e                          !a(size(a,1)-1,size(a,1)/2)=-7
 
      ! bloc 3
      mat(2,n+1) = b                         !a(2,(size(a,1)/2)+1) = 5
      mat(3,n+1) = -a                        !a(3,(size(a,1)/2)+1) = -4
      mat(n-1,n+1) = a                         ! a((size(a,1)/2)-1,(size(a,1)/2)+1) = 4
      mat(n,n+1) = -b                      !a(size(a,1)/2,(size(a,1)/2)+1) = -5
 
      mat(1,n+2) = -b                        !a(1,(size(a,1)/2)+2) = -5
      mat(3,n+2) = b                         !a(3,(size(a,1)/2)+2) = 5
      mat(4,n+2) = -a                        !a(4,(size(a,1)/2)+2) = -4
      mat(n,n+2) = a                       !a(size(a,1)/2,(size(a,1)/2)+2) = 4
 
      do i = n+3,2*n-2!(size(a,1)/2)-2
         mat(i-2-n,i) = a                  !a(i-2,(size(a,1)/2)+i) = 4
         mat(i-1-n,i) = -b                 !a(i-1,(size(a,1)/2)+i) = -5
         mat(i+1-n,i) = b                  ! a(i+1,(size(a,1)/2)+i) = 5
         mat(i+2-n,i) = -a                 !a(i+2,(size(a,1)/2)+i) = -4
      end do
 
      mat(1,2*n-1) = -a                      !a(1,size(a,1)-1) = -4
      mat(n-3,2*n-1) = a                   ! a((size(a,1)/2)-3,size(a,1)-1) = 4
      mat(n-2,2*n-1) = -b                    !a((size(a,1)/2)-2,size(a,1)-1) = -5
      mat(n,2*n-1) = b                     ! a(size(a,1)/2,size(a,1)-1) = 5
 
      mat(1,2*n) = b                       !a(1,size(a,1)) = 5
      mat(2,2*n) = -a                      ! a(2,size(a,1)) = -4
      mat(n-2,2*n) = a                   !a((size(a,1)/2)-2,size(a,1)) = 4
      mat(n-1,2*n) = -b                      !a((size(a,1)/2)-1,size(a,1)) = -5
 
 
 
      ! bloc 4
      mat(n+2,n+1) = -c                      !a((size(a,1)/2)+2,(size(a,1)/2)+1) = -3
      mat(2*n,n+1) = -c                    !a(size(a,1),(size(a,1)/2)+1) = -3
      do i = n+1,2*n!(size(a,1)/2)+1,2*size(a,1)
         mat(i,i) = d
      end do
      do i = n+2,2*n-1!(size(a,1)/2)+2,2*(size(a,1)/2)-1
         mat(i-1,i) = -c
         mat(i+1,i) = -c
      end do
      mat(n+1,2*n) = -c                    !a((size(a,1)/2)+1,size(a,1)) = -3
      mat(2*n-1,2*n) = -c                  !a(size(a,1)-1,size(a,1)) = -3 
      !print*,'mat = ',mat
    end subroutine matrice1
 
    subroutine transposemat(n,mat,tmat)
      use parametremod
      integer, intent(in) :: n
      real*8, dimension(2*n,2*n), intent(in) :: mat
      real*8, dimension(2*n,2*n), intent(out) :: tmat
      integer :: i,j
 
      !call matrice1(n,mat)
      do j=1,2*n
         do i=1,2*n
            tmat(j,i)=mat(i,j)
         end do
      end do
      !print*,tmat
    end subroutine transposemat
 
 
    subroutine matrice2(mat,n)
      use parametremod
      integer, intent(in) :: n
      real*8, dimension(2*n,2*n), intent(out) :: mat
      integer :: i
 
      ! Initialisation
      mat(:,:) = 0
      !print*,size(a,1)
 
      ! bloc 1
      do i = 1,n!size(a,1)/2
         mat(i,i) = 1
      end do
 
      !bloc 2
      mat(n+2,1)=-e                           !mat((size(a,1)/2)+2,1)=7
      mat(2*n,1)=e                        !a(size(a,1),1)=-7
      do i = 2,n-1!(size(a,1)/2)-1
         mat(i+n+1,i) = -e                    !a(i+(size(a,1)/2)+1,i) = 7
         mat(i+n-1,i) = e                     !a((size(a,1)/2)+i-1,i) = -7
      end do
      mat(n+1,n)=-e                         !a((size(a,1)/2)+1,size(a,1)/2)=7
      mat(2*n-1,n)=e                          !a(size(a,1)-1,size(a,1)/2)=-7
 
 
      ! bloc 3
      mat(2,n+1) = -b                         !a(2,(size(a,1)/2)+1) = 5
      mat(3,n+1) = a                        !a(3,(size(a,1)/2)+1) = -4
      mat(n-1,n+1) = -a                         ! a((size(a,1)/2)-1,(size(a,1)/2)+1) = 4
      mat(n,n+1) = b                      !a(size(a,1)/2,(size(a,1)/2)+1) = -5
 
      mat(1,n+2) = b                        !a(1,(size(a,1)/2)+2) = -5
      mat(3,n+2) = -b                         !a(3,(size(a,1)/2)+2) = 5
      mat(4,n+2) = a                        !a(4,(size(a,1)/2)+2) = -4
      mat(n,n+2) = -a                       !a(size(a,1)/2,(size(a,1)/2)+2) = 4
 
      do i = n+3,2*n-2!(size(a,1)/2)-2
         mat(i-2-n,i) = -a                  !a(i-2,(size(a,1)/2)+i) = 4
         mat(i-1-n,i) = b                 !a(i-1,(size(a,1)/2)+i) = -5
         mat(i+1-n,i) = -b                  ! a(i+1,(size(a,1)/2)+i) = 5
         mat(i+2-n,i) = a                 !a(i+2,(size(a,1)/2)+i) = -4
      end do
 
      mat(1,2*n-1) = a                      !a(1,size(a,1)-1) = -4
      mat(n-3,2*n-1) = -a                   ! a((size(a,1)/2)-3,size(a,1)-1) = 4
      mat(n-2,2*n-1) = b                    !a((size(a,1)/2)-2,size(a,1)-1) = -5
      mat(n,2*n-1) = -b                     ! a(size(a,1)/2,size(a,1)-1) = 5
 
      mat(1,2*n) = -b                       !a(1,size(a,1)) = 5
      mat(2,2*n) = a                      ! a(2,size(a,1)) = -4
      mat(n-2,2*n) = -a                   !a((size(a,1)/2)-2,size(a,1)) = 4
      mat(n-1,2*n) = b                      !a((size(a,1)/2)-1,size(a,1)) = -5
 
 
 
      ! bloc 4
      mat(n+2,n+1) = -c                      !a((size(a,1)/2)+2,(size(a,1)/2)+1) = -3
      mat(2*n,n+1) = -c                    !a(size(a,1),(size(a,1)/2)+1) = -3
      do i = n+1,2*n!(size(a,1)/2)+1,2*size(a,1)
         mat(i,i) = d
      end do
      do i = n+2,2*n-1!(size(a,1)/2)+2,2*(size(a,1)/2)-1
         mat(i-1,i) = -c
         mat(i+1,i) = -c
      end do
      mat(n+1,2*n) = -c                    !a((size(a,1)/2)+1,size(a,1)) = -3
      mat(2*n-1,2*n) = -c                  !a(size(a,1)-1,size(a,1)) = -3 
      !print*,'mat2 = ',mat
    end subroutine matrice2
  end module matricemod

Les coeff de la matrices

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
 
module parametremod
  implicit none
 
  real*8 :: a1,a2,a3,a,b1,b2,b,c,d1,d2,d,e,f,dx, dt, alpha, g,h
 
end module parametremod

les parametres:
-1.D0 1.D0 100 300 1.D0 1.2 9.8 0.01D0

[moomba]

Salut,

J'ai codé un GC il n'y a pas longtemps, pour résoudre l'équation de la chaleur. Si ca peut t'aider dans un premier temps, on regardera ton code ensuite si ca ne résoud pas ton problème. Note que la multiplication de la matrice se retrouve ici sous la forme :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
    ! Matrix Operation : q = A p
    do i=1,Nx(1)
       q(i) = (p(i+1)-2.0*p(i)+p(i-1)) * InvDxDx(1)
    end do

Le programme :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
program CG
 
 ! Conjugate Gradient
 !
 ! Benoît Leveugle, ASA
 !
 ! Solve 1D heat equation
 ! Solve A T = 0, with A the matrix corresponding to the laplacian
 !
 ! Heat equation : dT/dt = alpha d^2 T/dx^2
 ! When converged, dT/dt = 0, the equation becomes : 0 = alpha d^2 T/dx^2 <--> 0 = d^2 T/dx^2
 ! Numericaly, this means : (T(i+1)-2.0*T(i)+T(i-1)) * 1/(Dx*Dx) = 0 <--> A T = 0 with the matrix A, a diagonal symetric positiv definit matrix
 ! The matrix is like :
 ! -2/(dx*dx)	1/(dx*dx)	0		0		0		0		....
 ! 1/(dx*dx)	-2/(dx*dx)	1/(dx*dx)	0		0		0		....
 ! 0		1/(dx*dx)	-2/(dx*dx)	1/(dx*dx)	0		0		....
 ! 0		0		1/(dx*dx)	-2/(dx*dx)	1/(dx*dx)	0		....
 ! 0		0		0		1/(dx*dx)	-2/(dx*dx)	1/(dx*dx)	0
 ! etc
 
 implicit none
 
 integer, parameter :: ndim = 1
 real(8), dimension(1:ndim) :: Lx,Dx,InvDxDx
 integer, dimension(1:ndim) :: Nx
 real(8) :: beta, rho, rho0, gam
 integer :: n,i
 real(8), allocatable, dimension(:) :: T,Residut,p,q
 
 Nx(1) = 100
 
 Lx(:) = 0.1
 Dx(:) = Lx(:)/((Nx(:)-1)*1.0)
 InvDxDx(:) = 1.0/(Dx(:)*Dx(:))
 
 allocate(T(0:Nx(1)+1),Residut(0:Nx(1)+1),p(0:Nx(1)+1),q(0:Nx(1)+1))
 
 ! Initial Solution
 T(1:Nx(1)) = 0.0
 T(0) = 10.0
 T(Nx(1)+1) = -1.0
 
 ! Initialisation
 
 do i=1,Nx(1)
    Residut(i) = - (T(i+1)-2.0*T(i)+T(i-1)) * InvDxDx(1)
 end do
 
 p(:) = 0.0
 beta = 0.0
 
 rho = 0.0
 do i=1,Nx(1)
    rho = rho + Residut(i)*Residut(i)
 end do
 
 ! Iteration
 
 n = 0
 do
    n = n + 1
 
    do i=1,Nx(1)
       p(i) = Residut(i) + beta * p(i)
    end do
 
    ! Matrix Operation : q = A p
    do i=1,Nx(1)
       q(i) = (p(i+1)-2.0*p(i)+p(i-1)) * InvDxDx(1)
    end do
 
    gam = 0.0
    do i=1,Nx(1)
       gam = gam + p(i)*q(i)
    end do
    gam = rho / gam
 
    do i=1,Nx(1)
       T(i) = T(i) + gam * p(i)
    end do
 
    do i=1,Nx(1)
       Residut(i) = Residut(i) - gam * q(i)
    end do
 
    ! Test if residut is < 1.e-7 so that solution is converged
    if(maxval(abs(Residut(1:Nx(1)))) < 1e-7) exit ! Converged Solution
    write(*,*) "CG step ",n,"Residut",maxval(abs(Residut(1:Nx(1))))
 
    rho0 = rho ! Save old rho value
    rho = 0.0
    do i=1,Nx(1)
       rho = rho + Residut(i)*Residut(i)
    end do
 
    beta = rho / rho0
 
 
 end do
 
 write(*,*) "Solution converged in ", n ,"Iterations, with a res of ", maxval(abs(Residut(1:Nx(1))))
 
 ! Solution plot
 
 open(10,file="out-GC.dat")
 do i=1,Nx(1)
    write(10,*) Dx(1)*(i-1),T(i)
 end do
 close(10)
 
 deallocate(T,Residut,p,q)
 
end program CG

Mais bon, si ta matrice n'est pas définie positive, tu peut y aller franco et utiliser un BIGCSTAB, voir un BICTSGAB(2).

Cordialement