Buon pomeriggio. Ho un problema nel seguente codice. In particolare,nella funzione "iterazionejacobi",dovrebbero essere svolti un prodotto tra matrice e vettore e poi una somma di un altro vettore al risultato della moltiplicazione. Ora,la moltiplicazione è eseguita correttamente(sarebbe eseguita tra la matrice di iterazione C e il vettore x(k) al passo k-esimo. Esegue anche la successiva somma con il vettore(che è il vettore dei termini noti del sistema,ma cambiato,perché ogni elemento è stato diviso per il coeff. a
di ogni riga della matrice A iniziale,in modo da esplicitare la x nel sistema).Ma la esegue sbagliata,nel senso che è come se sbagliasse il vettore da sommare,nonostante il vettore d sia "preso" correttamente nel main. Non capisco il perché accada ciò. Eccovi il codice:
/*Programma C++ strutturato in funzioni che acquisisce da tastiera un numero intero 0<n<=15,gli elementi di una matrice A appartenente ad R^nxn,le componenti di un
vettore reale b appartenente ad R^n,e un numero reale 0<tol<=10^-4;successivamente approssima la soluzione del sistema lineare Ax=b con il metodo di Jacobi,
preventivando al massimo 100 iterazioni e utilizzando un criterio di arresto a posteriori che garantisca che la norma infinito dell'errore commesso all'arresto
della procedura sia sufficientemente piccola per garantire la precisione tol;infine stampa sul video il vettore "soluzione" e il numero di iterazioni compiute*/
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define N 15
#define MAXTOL 1.e-4
#define MAX_ITERATE 100
typedef double matrice[N][N];
typedef double vettore[N];
void letturaecontrollodimensione(int&);
void letturamatrice(matrice,int);
void stampamatrice(matrice,int);
void letturavettore(vettore,int);
void stampavettore(vettore,int);
double letturaprecisione(double&);
int algoritmojacobi(matrice,vettore,vettore,double,int);
void iterazionejacobi(matrice,vettore,vettore,vettore,int);
double costantelambda(matrice,int);
double modulo(double);
double normadelmassimo(vettore,int);
int main(){
printf("Programma C++ che legge un intero n tale che 0<n<=15,una matrice reale A di R^nxn,un vettore reale b di R^n e la tolleranza toll tale che 0<toll<=10^-4,approssima la \nsoluzione del sistema Ax=b con il metodo di Jacobi con max 100 iterazioni e un criterio di arresto a posteriori che garantisca che la norma infinito dell'errore all'arresto della procedura sia sufficientemente piccola da garantire la precisione toll,e stampa il vettore 'soluzione' e il numero di \niterazioni compiute\n");
int n;
matrice A;
vettore b;
letturaecontrollodimensione(n);
letturamatrice(A,n);
printf("La matrice dei coefficienti A[%d][%d] del sistema e':\n",n,n);
stampamatrice(A,n);
printf("\n");
letturavettore(b,n);
printf("Il vettore dei termini noti b[%d]del sistema e':\n",n);
stampavettore(b,n);
double lambda=costantelambda(A,n);
/*La costante di contrattività,lambda,deve essere tale che 0<=lambda<1(in modo tale che la matrice sia a diagonale
strettamente dominante)e perciò non può essere maggiore,o uguale,di (1-MAXTOL) (che è circa 1,essendo MAXTOL un numero molto
prossimo allo 0,al massimo 10^-4):*/
if(lambda>=1-MAXTOL) {printf("\nLa costante di contrattivita' lambda e'maggiore(o uguale) di 1:la matrice non e'a diagonale strettamente dominante.Termine del programma.\n");return 1;}
matrice C;/*matrice C di iterazione dell'algoritmo di Jacobi,ottenuta dalla matrice A esplicitando,per ogni riga i-esima di
A(cioè per ogni equazione del sistema),l'incognita x[i]*/
vettore d;/*vettore dei termini noti ottenuto da quello originario(b)dopo aver isolato l'incognita x[i]:ogni suo termine
sarà del tipo d[i]=b[i]/a[i][i]*/
//Costruzione della matrice di iterazione C e del vettore d:
/*L'elemento c[i][j] si ottiene spostando gli elementi della riga i al secondo membro dell'equazione (dove c'è il termine
noto),e dividendo entrambi i membri per il coefficiente A[i][i] del termine x[i](isolando quest'ultimo).Essendo un'equazione,
anche il termine i-esimo del vettore noto b viene diviso per A[i][i],ottenendo il termine i-esimo d[i] del vettore d:*/
for(int i=0;i<n;i++){
for(int j=0;j<n;j++){
int h=i;
if(i!=j) C[i][j]=-1*A[i][j]/A[h][h];
d[i]=b[i]/A[h][h];
}
}
printf("\nMatrice C:\n");
stampamatrice(C,n);
printf("\nVettore d:\n");
stampavettore(d,n);
double tol;
vettore xk; /*vettore di partenza(uguale al vettore nullo) per le iterazioni del metodo di Jacobi(il vettore al passo
successivo sarà il vettore b dei termini noti,essendo il vettore x(0) nullo)*/
for(int i=0;i<n;i++) xk[i]=1;
printf("\nVettore xk:\n");
stampavettore(xk,n);
double tau=letturaprecisione(tol);/*si legge la tolleranza tau alla quale dovrà arrivare il valore di
(lambda/(1-lambda))*||x(k-1)-x(k)||oo per far arrestare l'algoritmo di Jacobi*/
tau=(tau/lambda)*(1-lambda);/*quantità che sarà la tolleranza da non superare affinché sia possibile reiterare l'algoritmo
di Jacobi: in particolare,per il criterio di arresto a posteriori,l'errore al passo k dovrà essere minore
o uguale a (1-lambda)/lambda,moltiplicato per ||x(k-1)-x(k)||oo (cioè per la tolleranza tau)*/
int iterazioni=algoritmojacobi(C,d,xk,tau,n);//numero di iterazioni compiute nell'algoritmo di Jacobi
printf("\nIl vettore 'soluzione' x(k) che approssima la soluzione del sistema lineare Ax=b tramite l'algoritmo di Jacobi e':\n");
stampavettore(xk,n);
printf("\nIl numero di iterazioni compiute nell'applicazione dell'algoritmo di Jacobi al sistema Ax=b e'%d\n",iterazioni);
return 0;
}
void letturaecontrollodimensione(int& dim){
do{printf("Inserire la dimensione n del sistema lineare Ax=b,tale che 0<n<=15:\n");scanf("%d",&dim);} while(dim<0||dim>15);
return;
}
void letturamatrice(matrice a,int dim){
for(int i=0;i<dim;i++){
for(int j=0;j<dim;j++){
printf("Inserire l'elemento (reale)a[%d][%d] della matrice A:\n",i+1,j+1);scanf("%lf",&a[i][j]);}
}
return ;
}
void stampamatrice(matrice a,int dim){
for(int i=0;i<dim;i++){printf("\n");
for(int j=0;j<dim;j++){
printf("%lf\t",a[i][j]);}
}
return;
}
void letturavettore(vettore v,int dim){
for(int i=0;i<dim;i++)
{printf("Inserire l'elemento(reale) b[%d] del vettore b:\n",i+1);scanf("%lf",&v[i]);}
return;
}
void stampavettore(vettore v,int dim){
for(int i=0;i<dim;i++)
printf("%lf\t",v[i]);
return;
}
/*Lettura della precisione(tolleranza)per il massimo errore tra l'approssimazione della soluzione al passo k e quella al
passo k-1:*/
double letturaprecisione(double& t){
do{printf("\nInserire la tolleranza (reale)tol tale che 0<tol<=10^-4:\n");scanf("%lf",&t);} while(t<=0||t>MAXTOL);
return t;
}
//Funzione che implementa l'algoritmo di Jacobi(il vettore xk è quello delle soluzioni,ad ogni passo k-esimo):
int algoritmojacobi(matrice a,vettore v,vettore xk,double t,int dim){
vettore xkm1;//vettore delle soluzioni x,al passo (k-1)-esimo
int cont=0;/*contatore delle iterate dell'algoritmo di Jacobi(essendo questo "a posteriori",il numero di iterate compiute
è calcolato man mano che si va avanti con il procedimento)*/
do{ for(int j=0;j<dim;j++) xkm1[j]=xk[j]; //in x(k-1) si memorizza x(k),componente per componente
iterazionejacobi(a,v,xkm1,xk,dim); //Si fa l'iterazione di Jacobi(creando ogni vettore x(k) dal precedente,x(k-1) )
for(int j=0;j<dim;j++) xkm1[j]-=xk[j]; /*Al vettore x(k-1) (in cui è stato memorizzato x(k) ) si sottrae l'x(k) del
passo successivo: x(k) è quello che più si avvicina,man mano,alla soluzione
del sistema*/
//stampavettore(xk,dim);
printf("\t");
}
while((++cont<MAX_ITERATE)&&(normadelmassimo(xkm1,dim)>=t)); /*Si aumenta prima il cont(perché questo parte da 0)e poi si
controlla che sia minore di MAX_ITERATE(numero massimo di iterazioni voluto):questa condizione,unita alla limitazione che
la ||xkm1-x||oo (si è scritto,nel codice,xkm1,perché in questa è stata già memorizzata la differenza tra xkm1 e xk,dopo la
prima iterata)sia maggiore o uguale alla tolleranza t,costituisce il criterio di arresto a posteriori per l'algoritmo:si va
avanti cioè nel reiterare l'algoritmo di Jacobi finché valgono le due condizioni poste nel while*/
return cont;
}
/*Funzione che implementa l'iterazione di Jacobi,calcolando l'approssimazione xk della soluzione del sistema al passo
k-esimo,partendo da quella al passo precedente(il (k-1)-esimo),xmk1,che viene moltiplicata alla matrice di iterazione C,e
poi sommata al vettore d (C e d calcolati in precedenza dalla matrice A e dal vettore b dei termini noti):*/
void iterazionejacobi(matrice a,vettore v,vettore xkm1,vettore xk,int dim){
for(int i=0;i<dim;i++){
xk[i]=0;
for(int j=0;j<dim;j++)
xk[i]+=a[i][j]*xkm1[j]+v[i]; //x(k)=C*x(k-1)+d
}
printf("\nVettore alla iterazione k:\n");
stampavettore(xk,dim);
return;
}
/*Calcolo della costante di contrattività=lambda(max per h=1,.,dim della somma per j=1,..,dim(con j!=h)
di | a[h][j]/a[h][h] | ):*/
double costantelambda(matrice a,int dim){
double l;
double vettorelambda[dim];
l=vettorelambda[0];
for(int i=0;i<dim;i++){
for(int j=0;j<dim;j++){
if(j!=i) {vettorelambda[i]+=modulo(a[i][j]);}
vettorelambda[i]/=modulo(a[i][i]);
}
if(vettorelambda[i]>l) l=vettorelambda[i];
}
return l;
}
//Calcolo del modulo di un numero:
double modulo(double x){
if (x<0) x*=-1;
else x=x;
return x;
}
//Calcolo della norma infinito di un vettore(max per i=0,...,dim del modulo di w[i]):
double normadelmassimo(vettore w,int dim){
vettore assoluto;
double normainf=assoluto[0];
for(int i=0;i<dim;i++){
assoluto[i]=modulo(w[i]);
if(assoluto[i]>normainf) normainf=assoluto[i];
}
return normainf;
}