Salve a tutti, questo è il mio primo post, nel caso dovessi compiere qualche errore perdonatemi
Sto combattendo da qualche pomeriggio con un programma che deve trovare la soluzione di un sistema applicando Gauss e sostituzione all'indietro alla matrice associata A (b è il vettore dei termini noti), per poi calcolare la norma infinito del vettore residuo r=Ax-b.
Ecco il programma:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#define MAXN 100
int legginat(){
int n;
do{
printf("\n\nInserisci la dimensione della matrice\n\t");
scanf("%d",&n);
}while(n<0||n>=MAXN);
return(n);
}
void vett(int n,float b[]){
int i;
//vettore b dei termini noti
for(i=0;i<n;i++){
printf("\n\nInserisci la componente %d del vettore dei termini noti\n\t",i);
scanf("%f",&b[i]);
}
return;
}
void scambiocol(float a[MAXN][MAXN],int n,int i){
int j;
float app;
for(j=0;j<=n;j++){
app=a[j][i];
a[j][i]=a[j][i+1];
a[j][i+1]=app;
}
return;
}
void leggimat(int n, float a[MAXN][MAXN]){
int i,j;
//indici delle colonne
for(j=0;j<n;j++){
a[0][j]=j+1;
}
//coefficienti della matrice A
for(i=1;i<n+1;i++){
for(j=0;j<n;j++){
printf("\n\nInserisci la componente %d,%d della matrice\n\t",i-1,j);
scanf("%f",&a[i][j]);
}
}
//stampa della matrice inserita
printf("\n\nMatrice iniziale:\n");
for(i=0;i<=n;i++){
for(j=0;j<n;j++){
printf("\t%f",a[i][j]);
}
printf("\n");
}
return;
}
int massimo(float a[MAXN][MAXN],int n,int k){
int i,imax;
float app;
app=a[k+1][k];
imax=k+1;
for(i=k+1;i<=n;i++) {
if(fabs(a[i][k])>app)
imax=i;
}
return(i);
}
void scambiorighe(float a[MAXN][MAXN],int n,int imax,int i){
int j;
float app;
for(j=0;j<n;j++){
app=a[i][j];
a[i][j]=a[imax][j];
a[imax][j]=app;
}
return;
}
void gauss(int n,float a[MAXN][MAXN],float b[]){
int i,j,k;
float c,imax,m;
for(k=0;k<n-1;k++){
imax=massimo(a,n,k);
if(imax=!k+1) scambiorighe(a,n,imax,k+1);
for(i=k+2;i<=n;i++){
//controllo pivot
if(a[k+1][k]==0) scambiocol(a,n,k);
m=a[i][k]/a[k+1][k]; printf("\n\nValore %d di m: %f",k,m); //stampa di controllo
for(j=k;j<n;j++){
a[i][j]=a[i][j]-m*a[k+1][j];
}
b[i-1]=b[i-1]-m*b[k];
}
//stampa della matrice parziale
printf("\n\nMatrice parziale, passaggio %d:\n",k);
for(i=0;i<=n;i++){
for(j=0;j<n;j++){
printf("\t%f",a[i][j]);
}
printf("\n");
}
}
//stampa di U
printf("\n\nMatrice triangolarizzata:\n");
for(i=0;i<=n;i++){
for(j=0;j<n;j++){
printf("\t%f",a[i][j]);
}
printf("\n");
}
//singolarità
for(i=0;i<n;i++){
if(a[i+1][i]==0){
printf("\n\nLa matrice e' SINGOLARE, non esiste un'unica soluzione");
leggimat(n,a);
}
}
return;
}
void sost_indietro(int n,float U[MAXN][MAXN], float b[MAXN],float x[MAXN]){
int i,j,k;
float som;
x[n-1]=b[n-1]/U[n][n-1];
for(k=n-1;k>0;k--){
som=0;
for(j=k+1;j<n;j++){
som+=(U[k][j]*x[j]);
}
x[k]=(1./U[k][k])*(b[k]-som);
}
return;
}
void stampa(int n,float x[]){
printf("\n\nLa soluzione del sistema associato alla matrice e':\n");
for(int i=0;i<n;i++){
printf("%f\n",x[i]);}
return;
}
void prods(int n,float a[MAXN][MAXN],float x[],float prod[]) {
int i,j;
//inizializzazione
for(i=0;i<n;i++){
prod[i]=0;
}
//prodotto righe per colonne Ax
for(j=0;j<=n;j++){
for(i=0;i<n;i++){
prod[i]=prod[i]+(a[j][i]*x[i]);
}
}
return;
}
float norma(int n,float a[MAXN][MAXN],float b[],float x[],float r[]){
float prod[MAXN],res=0;
int i;
for(i=0;i<n;i++){
prods(n,a,x,prod);
r[i]=b[i]-prod[i];
}
for(i=0;i<n;i++){
if (res<r[i])
res=r[i];
}
return(res);
}
int main (){
int n;
float a[MAXN][MAXN],b[MAXN],x[MAXN],r[MAXN],res;
n=legginat();
leggimat(n,a);
vett(n,b);
gauss(n,a,b);
sost_indietro(n,a,b,x);
stampa(n,x);
res=norma(n,a,b,x,r);
printf("\n\n La norma infinito e' %f",res);
return (1);
}
Ho provato il programma con
A:
12 1 2 3
1 12 2 3
1 2 12 3
1 2 3 12
b: 18, 18, 18, 18
In teoria la soluzione dovrebbe essere un vettore x= 1,1,1,1
Ecco l'output:
Valore 0 di m:12
valore 0 di m:1
Valore 0 di m:1
Matrice parziale, passaggio 0:
1 2 3 4
1 12 2 3
0 -143 -22 -33
0 -10 10 0
0 -10 1 9
Valore 1 di m:-0.833333
Valore 1 di m:-0.833333
Matrice parziale, passaggio 1:
1 2 3 4
0 -143 -22 -33
1 12 2 3
0 -0 11.666667 2.5
0 -0 2.666667 11.5
Valore 2 di m: -0.121212
Matrice parziale, passaggio 2:
1 2 3 4
0 -0 11.666667 2.5
1 12 2 3
0 -143 -22 -33
0 -0 -0 7.5
Matrice triangolarizzata:
1 2 3 4
0 -0 11.666667 2.5
1 12 2 3
0 -143 -22 -33
0 -0 -0 7.5
La matrice è SINGOLARE, non esiste un'unica soluzione
Inserisci la componente 0,0 della matrice
Probabilmente c'è qualche problema serio nella funzione gauss, ma non riesco a trovarlo.. in più proprio non riesco a capire come mai il valore di m venga calcolato più di una volta!
HELP