Rundungsfehler?



  • Hallo alle miteinander.

    Ich hoffe ich bin im richtigen Forum.
    Ich habe die Aufgabe die Maxwellschen Gleichungen mit einem expliziten Differenzenverfahren zu lösen. Das habe ichmal in C++ gemacht und wollte fragen ob jemand weiß warum mein ergebnis so komisch ist.
    Nach 10 Zeitschritten sieht meine Lösung so aus (kann man eigentlich auch Bilder direkt ins Forum einbinden?):

    http://www.freeimagehosting.net/image.php?812ee5249e.png

    Das ist noch so wie ich es gehofft habe. Nach 20 zeitschritten sieht das ganze aber nicht mehr so schön Sinus-förmig aus:

    http://www.freeimagehosting.net/image.php?c0142d5bef.png

    Kann sich das jemand erklären? Warum ist die Lösung in der Nähe der Null viel kleiner und warum artet das nach außen hin so aus? Ich habe nicht wirklich viel Ahnung von Computerarithmetik aber ich vermute dass das Rundungsfehler sind, die am Rand größer sind, weil dort 2 benachbarte Fließkommazahlen da weiter auseinander liegen.
    Meine Daten haben den Datentyp long double.

    Oder hättet ihr eine Begründung?
    Den Code poste ich bewusst nicht mit, weil das den Rahmen sprengen würde.
    Meine zweite Erklärung wäre, dass das Verfahren sehr instabil ist und sich Fehler schnell aufsummieren könnten. Aber warum dann nur so am Rand?

    Grüße

    Max



  • Max3000 schrieb:

    Ich habe die Aufgabe die Maxwellschen Gleichungen mit einem expliziten Differenzenverfahren zu lösen.

    Das sagt mir beides nichts. Was für ein mathematisches Problem ist das und was meinst Du mit Differenzenverfahren? Willst Du ein Anfangswertproblem (Differentialgleichungen) lösen?

    Vielleicht fragst Du lieber in einem Mathe/Numerik-spezifischen Forum nach...



  • Du solltest das Problem wirklich erstmal beschreiben. Maxwell-Gleichungen, schön und gut... wie sind die Rand-/Anfangsbedingungen? Formuliere mal in Worten, was Du implementiert hast, wenn Du den Code schon nicht zeigen magst. Hast Du die Stabilität Deiner Methode gecheckt? Nach Rundungsfehlern sieht das nämlich nicht aus.



  • Ok also die Gleichungen lauten:

    μH_xt=E_zy\mu\frac{\partial H\_x}{\partial t}=-\frac{\partial E\_z}{\partial y}

    μH_yt=E_zx\mu\frac{\partial H\_y}{\partial t}=\frac{\partial E\_z}{\partial x}

    ϵE_zt=H_yxHxy\epsilon\frac{\partial E\_z}{\partial t}=\frac{\partial H\_y}{\partial x}-\frac{\partial H_x}{\partial y}

    , mit den Anfangsbedingungen

    $H_x^0(x,y)=0$\newline $H_y^0(x,y)=-sin(2\pi x/500)*(\epsilon/\mu)^{1/2}$\newline $E_z^0(x,y)=sin(2\pi x/500)$\newline

    Am Rand liegen periodische Randbedingungen vor.

    Das ganze sollten wir mit dem sogenannten Lax-Friedrich-Verfahren lösen, was den neuen Zeitschritt aus dem letzten wie folgt berechnet:

    $u_{i,j}^{n+1}=1/4(u_{i+1,j}^n+u_{i-1,j}^n+u_{i,j+1}^n+u_{i,j-1}^n)$\newline $-\frac{\Delta t}{2*\Delta x}*(F_{x,i+1,j}^n-F_{x,i-1,j}^n+F_{y,i,j+1}^n-F_{y,i,j-1}^n)$

    ,wobei Fx und Fy irgendwie definiert sind, aber das ist ja nur sinnloses einsetzen.

    Und zwar ein bisschen Code noch:

    // Funktoren für die Anfangsbedingung:
    struct initial_value_Hx
    {
      long double operator()(long double x=0, long double y=0) 
      {
        return 0;
      }
    };
    
    struct initial_value_Hy
    {
      // Define static constants
      const long double mu;
      const long double eps;
    
      initial_value_Hy() : mu(4*M_PI*1e-7), eps(8.842e-12) {}
    
      long double operator()(long double x=0, long double y=0) 
      {
        return sin(2.*M_PI/500.*x)*(-sqrt(eps/mu));
      }
    };
    
    struct initial_value_Ez
    {
      long double operator()(long double x=0, long double y=0) 
      {
        return sin(2.*M_PI/500.*x);
      }
    };
    
    // Initialisierung der Anfangslösung mit eben diesen Funktoren
      template <typename HxFunctor, typename HyFunctor, typename EzFunctor>
      Solution(HxFunctor &hx, HyFunctor &hy, EzFunctor &ez) 
        : Hx(N, N), Hy(N, N), Ez(N, N), Dx(2000/N)
      {
        for(unsigned i=0; i<N; ++i)
          for(unsigned j=0; j<N; ++j)
    	{
    	  Hx(i, j) = hx(get_x_coord(i, N), get_y_coord(j, N));
    	  Hy(i, j) = hy(get_x_coord(i, N), get_y_coord(j, N));
    	  Ez(i, j) = ez(get_x_coord(i, N), get_y_coord(j, N));
    	}
      }
    
    // Und hier die Berechnung des neuen aus dem alten Zeitschritt. Der Rand ist periodisch und muss extra behandelt werden
    
      void fromPreviousTimestep(Solution &prev) 
      {
        // Calculate interior points
        for(unsigned i=1; i<N; ++i)
          for(unsigned j=1; j<N; ++j)
    	{
    	  Hx(i, j) = 1/4*(prev.Hx(i-1,j)+prev.Hx(i,j-1)+prev.Hx((i+1)%N,j)+prev.Hx(i,(j+1)%N))
    	    + Dt/(2*Dx)*(prev.Fx_0(i-1,j) + prev.Fy_0(i,j-1) - (prev.Fx_0((i+1)%N,j) + prev.Fy_0(i,(j+1)%N)));
    	  Hy(i, j) = 1/4*(prev.Hy(i-1,j)+prev.Hy(i,j-1)+prev.Hy((i+1)%N,j)+prev.Hy(i,(j+1)%N))
    	    + Dt/(2*Dx)*(prev.Fx_1(i-1,j) + prev.Fy_1(i,j-1) - (prev.Fx_1((i+1)%N,j) + prev.Fy_1(i,(j+1)%N)));
    	  Ez(i, j) = 1/4*(prev.Ez(i-1,j)+prev.Ez(i,j-1)+prev.Ez((i+1)%N,j)+prev.Ez(i,(j+1)%N))
    	    + Dt/(2*Dx)*(prev.Fx_2(i-1,j) + prev.Fy_2(i,j-1) - (prev.Fx_2((i+1)%N,j) + prev.Fy_2(i,(j+1)%N)));
    	}
    
        // Calculate Boundary points
        Hx(0,0) = 1/4*(prev.Hx(1,0)+prev.Hx(N-1,0)+prev.Hx(0,1)+prev.Hx(0,N-1)) 
          + Dt/(2*Dx)*(prev.Fx_0(N-1,0)+prev.Fy_0(0,N-1)-(prev.Fx_0(1,0)+prev.Fy_0(0,1)));
        Hy(0,0) = 1/4*(prev.Hy(1,0)+prev.Hy(N-1,0)+prev.Hy(0,1)+prev.Hy(0,N-1)) 
          + Dt/(2*Dx)*(prev.Fx_1(N-1,0)+prev.Fy_1(0,N-1)-(prev.Fx_1(1,0)+prev.Fy_1(0,1)));
        Ez(0,0) = 1/4*(prev.Ez(1,0)+prev.Ez(N-1,0)+prev.Ez(0,1)+prev.Ez(0,N-1)) 
          + Dt/(2*Dx)*(prev.Fx_2(N-1,0)+prev.Fy_2(0,N-1)-(prev.Fx_2(1,0)+prev.Fy_2(0,1)));
    
        for(unsigned i=1; i<N; ++i) {
          Hx(0,i) = 1/4*(prev.Hx(0,(i+1)%N)+prev.Hx(0,i-1)+prev.Hx(1,i)+prev.Hx(N-1,i)) 
    	+ Dt/(2*Dx)*(prev.Fx_0(N-1,i)+prev.Fy_0(0,i-1)-(prev.Fx_0(1,i)+prev.Fy_0(0,(i+1)%N)));
          Hy(0,i) = 1/4*(prev.Hy(0,(i+1)%N)+prev.Hy(0,i-1)+prev.Hy(1,i)+prev.Hy(N-1,i)) 
    	+ Dt/(2*Dx)*(prev.Fx_1(N-1,i)+prev.Fy_1(0,i-1)-(prev.Fx_1(1,i)+prev.Fy_1(0,(i+1)%N)));
          Ez(0,i) = 1/4*(prev.Ez(0,(i+1)%N)+prev.Ez(0,i-1)+prev.Ez(1,i)+prev.Ez(N-1,i)) 
    	+ Dt/(2*Dx)*(prev.Fx_2(N-1,i)+prev.Fy_2(0,i-1)-(prev.Fx_2(1,i)+prev.Fy_2(0,(i+1)%N)));
        }
        for(unsigned i=1; i<N; ++i) {
          Hx(i,0) = 1/4*(prev.Hx((i+1)%N,0)+prev.Hx(i-1,0)+prev.Hx(i,1)+prev.Hx(i,N-1)) 
    	+ Dt/(2*Dx)*(prev.Fx_0(i-1,0)+prev.Fy_0(i,N-1)-(prev.Fx_0((i+1)%N,0)+prev.Fy_0(i,1)));
          Hy(i,0) = 1/4*(prev.Hy((i+1)%N,0)+prev.Hy(i-1,0)+prev.Hy(i,1)+prev.Hy(i,N-1)) 
    	+ Dt/(2*Dx)*(prev.Fx_1(i-1,0)+prev.Fy_1(i,N-1)-(prev.Fx_1((i+1)%N,0)+prev.Fy_1(i,1)));
          Ez(i,0) = 1/4*(prev.Ez((i+1)%N,0)+prev.Ez(i-1,0)+prev.Ez(i,1)+prev.Ez(i,N-1)) 
    	+ Dt/(2*Dx)*(prev.Fx_2(i-1,0)+prev.Fy_2(i,N-1)-(prev.Fx_2((i+1)%N,0)+prev.Fy_2(i,1)));
        }
    
      }
    


  • Probier es doch vielleicht mal mit einer trivialen Problemstellung aus, z.B. H=E=0. Dann dürfte eigentlich nichts passieren und man kann ausschließen, dass Du evtl. irgendwo über den Speicher hinaus liest o.ä. Danach würde ich mal die Flüsse checken und vor allem erstmal schauen, was passiert, wenn man den Zeitschritt deutlich verkleinert.



  • Also an dem $$\Delta x$ und Δt\Delta t$ hab ich schon viel rumgespielt und ab einem gewissen Schritt artet das immer so aus.
    Aber das E=H_x=0 probier ich mal aus. Das ist eine super Idee. Danke.



  • Falscher Datentyp... Ich tippe drauf das die Abweichungen durch die Ungenauigkeit des Datentyp "double" entstehen. Hierdurch werden ständig Werte leicht abgerundet in der Berechnung. Double kann zwar einen sehr großen Zahlenbereich abbilden, macht dies aber indem nur eine Teilmenge dieses Bereiches abgebildet wird. Alle nicht abbildbaren Zahlen werden automatisch auf die nächste abbildbare abgerundet.

    Im Prinzip das gleiche Konzept wie bei Integer, wo ja praktisch alle Nachkomma-Stellen verschluckt werden, da ein double aus zwei Integer-Werten aufgebaut wird, nämlich Mantisse und Exponent.


Log in to reply