Momentaneamente sto lavorando al metodo di risoluzione dei sistemi lineari tramite Jacobi.
Ho scritto il codice in C ma non riesco a farlo funzionare. Tenete conto che sono alle primissime armi con il linguaggio C. Vi posto il codice, se qualcuno ha suggerimenti sono ben accetti!
Il compilatore comunque non da errori di sintassi!
- Codice:
//Algoritmo di Jacobi [Chiamata MatLab: [x,its]=jacobi(A,b,x_0,maxit,tol)]
/* L'algoritmo deve prendere in input una matrice quadrata A con elementi della diagonale non nulli,
un vettore (termine noto) b, una approssimazione della soluzione iniziale x_0, e i parametri d'arresto
dati da maxit= numero massimo di iterezioni consentite e tol=tolleranza errore.
L'algoritmo deve restituire il vettore soluzione x e il parametro its delle iterazioni compiute.*/
#include<stdio.h>
#include<math.h>
#define SIZE 50
double molt(size_t dim, double mat[dim][dim],double vet[dim]);
double norma(size_t dim,double vet[dim]);
double molt_mat(size_t dim,double mat1[dim][dim],double mat2[dim][dim]);
double somma_vet(size_t dim, double vec1[dim],double vec2[dim]);
int main(){
size_t dim;
int i,j,its,maxit;
double tol,norm;
printf("-----METODO JACOBI-----\n");
printf("Si inseriscano i dati nella maniera richiesta.\n");
printf("Si inserisca la dimensione della matrice A.\n");
scanf("%d", &dim);
double A[dim][dim];
double b[dim];
double x0[dim];
double x[dim];
double err[dim];
printf("Si inseriscano gli elementi della matrice A nell'ordine rischiesto.\n");
for (i=0; i<dim; i++){
for(j=0;j<dim;j++){
printf("Inserisci l'elemento di posto %d , %d : \n",i+1,j+1);
scanf("%lf", &A[i][j]);
}
}
//controllo sugli elementi di A che nn ci siano elementi della diag nulli
for (i=0;i<dim;i++){
if(A[i][i]==0.00){
printf("Errore: non e' possibile applicare il metodo di Jacobi.");
return 0;
}
}
printf("Si inseriscano gli elementi del termine noto b nell'ordine richiesto.\n");
for (i=0; i<dim; i++){
printf("Inserisci l'elemento di posto %d: \n",i+1);
scanf("%lf", &b[i]);
}
printf("Si inseriscano gli elementi del vettore iniziale x_0.\n");
for (i=0; i<dim; i++){
printf("Inserisci l'elemento di posto %d: \n",i+1);
scanf("%lf", &x0[i]);
}
printf("Si inserisca il numero di iterazioni massime.\n");
scanf("%d",&maxit);
printf("Si inserisca la tolleranza.\n");
scanf("%g", &tol);
printf("toll=%lg",&tol);
//---------------------------------------------
//creo i pezzi per l'iterazione
//inizializzo x a x0
for(i=0;i<dim;i++){
x[i]=x0[i];
}
//calcolo la norma iniziale
for(i=0;i<dim;i++){
err[i]=b[i]-molt(dim,A,x);
}
norm=norma(dim,err);
printf("norma %g",norm);
//creo matrice diagonale
double D[dim][dim];
for(i=0;i<dim;i++){
for(j=0;j<dim;j++){
if(i==j){
D[i][j]=A[i][i];
}
else{
D[i][j]=0.00;
}
}
}
//Creo D^-1
double D_inv[dim][dim];
for (i=0;i<dim;i++){
D_inv[i][i]=1/D[i][i];
}
//creo N=A-D
double N[dim][dim];
for (i=0;i<dim;i++){
for(j=0;j<dim;j++){
if(i==j){
N[i][j]=0.00;
}
else{
N[i][j]=A[i][j];
}
}
}
//creo T=D^-1*N
double T[dim][dim];
T[dim][dim]=molt_mat(dim,D_inv,N);
//creo S=D_inv*b
double S[dim];
S[dim]=molt(dim,D_inv,b);
//ITERAZIONE JACOBI
int iter=0;
while(iter<maxit && norm>tol){
double temp[dim];
temp[dim]=molt(dim,T,x);
x[dim]=somma_vet(dim,temp,S);
iter++;
for(i=0;i<dim;i++){
err[i]=b[i]-molt(dim,A,x);
}
norm=norma(dim,err);
}
//presentazione risultati
printf("La soluzione calcolata e':\n");
for (i=0;i<dim;i++){
printf("%f ",x[i]);
printf("\n");
}
printf("\n");
printf("le iterazioni fatte sono %d \n",iter);
return 0;
}
//_________FUNZIONI A SUPPORTO___________
double molt(size_t dim, double mat[dim][dim], double vet[dim]){
int i,j;
double sum;
double RIS[dim];
for(i=0;i<dim;i++){
sum=0.00;
for(j=0;j<dim;j++){
sum+=mat[i][j]*vet[j];
}
RIS[i]=sum;
}
return RIS[dim];
}
//____________
double norma(size_t dim, double vet[dim]){
int i;
double sum,norm;
sum=0.00;
for(i=0;i<dim;i++){
sum+=pow(vet[i],2);
}
norm=fabs(sqrt(sum));
return norm;
}
//__________
double molt_mat(size_t dim,double mat1[dim][dim],double mat2[dim][dim]){
double RIS[dim][dim];
double sum;
int n,i,j,k;
for(i=0;i<dim;i++){
for(j=0;j<dim;j++){
sum=0.00;
for(n=0;n<dim;n++){
sum+=mat1[i][n]*mat2[n][j];
}
RIS[i][j]=sum;
}
}
return RIS[dim][dim];
}
//___________
double somma_vet(size_t dim, double vec1[dim],double vec2[dim]){
double RIS[dim];
int i;
for(i=0;i<dim;i++){
RIS[i]=vec1[i]+vec2[i];
}
return RIS[dim];
}
