Ciclo for o non ciclo for: questo è il problema!

[fatallegra]

Buonasera a tutti!
Vi pongo un problema in cui mi sono imbattuta nei tre giorni appena trascorsi: sto lavorando ad un programma scritto in C++ mediante il quale ricavare il raggio spettrale di una matrice al variare di un certo coefficiente "tau" (in particolare devo verificare la stabilità del metodo FTCS) e mi sono accorta che ottengo risultati diversi se faccio variare "tau" con un ciclo for o se lo faccio variare "manualmente", cioè cambiando di volta in volta il suo valore prima di compilare ogni volta il programma.
Chiedo scusa a quanti stanno leggendo il mio post, come avrete capito sono una principiante, quindi mi scuso soprattutto per il linguaggio che sto usando, probabilmente inappropriato .
Il programma è il seguente (non so se possa essere utile):

#include
#include
#include
#include

using namespace std;

int main()
{
int N=101;
double I[101][101], D[101][101], A[101][101], B[101][101];
double sum;
double rho;
double tau=0.000000001;
double h=0.01;
double t=(h*h)/2;


ofstream file("B.dat");

for(int i=0; i {
for(int j=0; j {
if(i==j)
{
I[j]=1;
}
else
{
I[j]=0;
}
}
}

for(int i=0; i {
for(int j=0; j {
if(j==i+1)
{
if(j==0)
{
D[j]=0;
}
else
{
D[j]=1;
}
}
else if(j==i)
{
if(j==0 || j==N-1)
{
D[j]=0;
}
else
{
D[j]=-2;
}
}
else if(j==i-1)
{
if(j==0)
{
D[j]=0;
}
else
{
D[j]=1;
}
}
else
{
D[j]=0;
}
}
}




while(tau<=100)
{

for(int i=0; i {
for(int j=0; j {
if(i==j)
{
double c=tau/(2*t);
A[j]=I[j]+c*D[j];
B[j]=A[j]*A[j];
sum+=B[j];
}
}
}
rho=sqrt(sum);
file << tau << "\t" << rho << endl;
tau*=10;
}

fflush(stdin);
return 0;
}

NON utilizzando il ciclo for ottengo valori del raggio spettrale che si aggirano tutti intorno ai 10, finchè non si ha -giustamente- divergenza, mentre se USO il ciclo for l'ordine di grandezza è lo stesso, però i valori sono diversi tra loro.
Il mio quesito è: qualcuno potrebbe avvertirmi se sto sbagliando qualcosa e se sì per quale motivo?
Vi ringrazio anticipatamente,

Francesca

[claudiocarcaci]

To consiglio di usare il tag

code

per rendere più leggibile il codice.

Comunque, intendi dire che se scrivi qualcosa tipo:

for(...)
   for(...)
      calcolo sum
calcolo rho
incremento tau

for(...)
   for(...)
      calcolo sum
calcolo rho
incremento tau

per 10 volte consecutive anzichè usare il ciclo while ottieni $ tau = 10 $ ?
Invece col ciclo while ottieni valori compresi tra -10 e +10?

[vict85]

Ho cambiato qualcosa e commentato. Dando principalmente qualche suggerimento per semplificare il codice. L'output è esattamente lo stesso del tuo.

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <fstream>

using namespace std;

int main()
{
    unsigned int const N = 101; // meglio metterlo in const.
    //hai definito una variabile per la definizione, allora usala. Inoltre ho mandato I e D a zero.
    double I[N][N] = { 0 }, D[N][N] = { 0 }, A[N][N], B[N][N];
    //somma non inizializzata. ATTENZIONE la prossima volta. Non sono sicuro però vada inizializzato qui.
    double sum=0.0;
    double rho;
    double tau = 1E-9; // scrittura equivalente ma penso più espressiva.
    double const h = 0.01; // penso sia anche questa costante.
    double t = (h*h)/2;

    // Per definire su I la matrice identità ho cambiato il tuo codice mettendo tutto a zero
    // nella definizione e ora modificando solo gli elementi sulla diagonale.
    // Molto più semplice, non ti sembra?
    for(unsigned int i=0; i!=N; ++i)
    {
      I[i][i] = 1.0;
    }

    // La definizione di D è difficile da capire così e troppo lunga.
    // Ho messo a 0 tutto e ora scrivo sopra solo le cose essenziali.
    // Alcune tue condizioni tra l'altro hanno poco senso come chiedersi se j==0 sapendo che j == i+1>= 1
    // Il mio codice dà per scontato che N >= 3. Non penso sia particolarmente restrittivo.
    // In caso bisognerebbe aggiungere un test su N (che verrebbe selezionato in fase di compilazione).
    // Controlla se ho preso in considerazione tutti i casi. Tieni conto che tutto il resto è 0.
    D[0][1] = 1.0;
    D[1][1] = -2.0;
    D[1][2] = 1.0;
    for(unsigned int i=2; i!=N-1; ++i)
    {
      D[i][i-1] = 1.0;
      D[i][i] = -2.0;
      D[i][i+1] = 1.0;
    }
    D[N-1][N-2] = 1.0;

    // meglio aprire i file dove li usi. Tanto non ha alcun limite a dove definire una variabile.
    // Eventualmente potevi creare prima ed aprire dopo.
    ofstream file("B.dat");

    while(tau <= 100)
    {
      // questa parte non ha senso. Stai facendo una operazione solo nel caso in cui sono uguali,
      // a quel punto creare 4 matrici era inutile... Penso che in questo ciclo ci sia un errore grave.
      for(unsigned int i=0; i!=N; ++i)
      {
        for(unsigned int j=0; j!=N; ++j)
        {
          if(i==j)
          {
            double const c=tau/(2*t);
            A[i][j]=I[i][j]+c*D[i][j];
            B[i][j]=A[i][j]*A[i][j];
            sum+=B[i][j];
          }
        }
      }


      rho=sqrt(sum);
      file << tau << "\t" << rho << endl;
      tau*=10.0;
    }

    // è inutile usare fflush su stdin anche considerando che non usi stdio ma iostream (non so come fa a compilare).
    // il comportamento di fflush su uno stream di input non è neanche uniforme a seconda della libreria che usi
    // non essendo specificato nello standard. Chi te l'ha insegnato come minimo dovrebbe ripassarsi lo standard C.
    // Attenzione inoltre che non hai chiuso il file.
    file.close();

    return 0;
}

Come puoi leggere dal codice che ti ho postato c'é una parte che non ho modificato (tranne un paio di preferenze stilistiche mie) e che a mio avviso è sbagliata:

for(unsigned int i=0; i!=N; ++i)
{
  for(unsigned int j=0; j!=N; ++j)
  {
    if(i==j)
    {
      double const c=tau/(2*t);
      A[i][j]=I[i][j]+c*D[i][j];
      B[i][j]=A[i][j]*A[i][j];
      sum+=B[i][j];
    }
  }
}

Tieni conto che quel codice è equivalente a

for(unsigned int i=0; i!=N; ++i)
{
  double const c=tau/(2*t);
  A[i][i]=I[i][i]+c*D[i][i]; // o anche = 1+c*D[i][i]; fare un if è inutile... Il più delle volte è 1-2c
  B[i][i]=A[i][i]*A[i][i];
  sum+=B[i][i];
}

Immagino invece che il tutto vada espresso in forma matriciale come A = cD + I, B=A*A (component-wise o come matrice) e sum come la somma degli elementi di B... Se è così il calcolo è ben diverso.

Tieni conto che esistono modi comodi e compatti per salvare matrici tridiagonali come sono I, D e immagino A e B. Inoltre non è detto che convenga avere A invece che semplicemente sommare per A*A... Il tutto, tra l'altro con più cicli più comodi come quello che ho usato per scrivere D.

Comunque il tuo problema con il ciclo è senza dubbio il problema di inizializzazione di sum... Ma secondo me il programma è probabilmente sbagliato in ogni caso.