Entwicklung von Differential-Gleichungen simulieren

[ghostwhisperer]

Hallo!
Ich hab mal ein grundlegend mathematisches Problem:

Wie simuliert man eine Differential-Gleichung? Sprich, wie funktioniert eine numerische Entwicklung?
Ganz einfaches Modell: eine linkslaufende (+k) und eine rechtslaufende Welle (-k)
sollten sich zu einer Stehwelle summieren.

Ich hab jetzt das Problem, dass die numerische zeitliche Entwicklung asymmetrisch ist. Und zwar egal ob in Excel oder in SciLab. Es liegt also am mathematisch-numerischen Ansatz. Entweder beide Wellen evolvieren in plus oder in minus-Richtung. Daran ist im folgenden Pseudocode nur ein Vorzeichen schuld. Ob ich erste und zweite Ableitung auf die "Integrale" addiere oder subtrahiere. Im Grunde sehen die Wellen in beiden Fällen ok aus, bis auf diese zeitliche Asymmetrie..
Der Pseudocode in SciLab (abgesehen von initialisieren und Schleifen usw):

DAs ganze ist nur eindimensional. Was mache ich falsch beim "integrieren"?
- das ganze auf einer "geschlossenen Kurve". Die Ränder sind verknüpft.

fi: erste Ableitung
ffi: zweite Ableitung
v: erste Summation "Geschwindigkeit"
s: zweite Summation "Elongation"


for i=1:nFrames

for j=2:1999
fi(j)=(u(j+1)-u(j-1))/2/0.006;
end
fi(1)=(u(2)-u(2000))/2/0.006;
fi(2000)=(u(1)-u(1999))/2/0.006;

for j=2:1999
ffi(j)=(fi(j+1)-fi(j-1))/2/0.006;
end
ffi(1)=(fi(2)-fi(2000))/2/0.006;
ffi(2000)=(fi(1)-fi(1999))/2/0.006;

for j=1:2000
v(j) = ffi(j)*0.006+fi(j);
end

for j=1:2000
s(j)= v(j)*0.006+s(j);
u(j)=s(j);
end
end


Hab bisher keine Arbeit in der Richtung gefunden.. DANKE!

GRüße, ghosti

[Geku]

Ich würde die Simulation in C wie folgt lösen:

PHP-Code:

#include <stdio.h>
#include <math.h>

#define N 100 // Anzahl der Diskretisierungspunkte im Raum
#define M 1000 // Anzahl der Diskretisierungspunkte in der Zeit
#define A 1.0 // Amplitude der Welle
#define k 1.0 // Wellenzahl
#define v 1.0 // Geschwindigkeit der Welle


int main()
{
    double y[N][M]; // Amplitude an den Diskretisierungspunkten (Raum/Zeitpunkte)
    double dx = 1.0 / N; // räumliche Schrittweite
    double dt = 0.001; // zeitliche Schrittweite
    double c = v * dt / dx; // CFL-Zahl
    double t, x;

    // Initialisierung der Amplitude am Anfang der Simulation

    for (int i = 0; i < N; i++) 
    {
        y[i][0] = A * cos(k * i * dx);
    }

    // Berechnung der Amplitude zu jedem Zeitpunkt

    for (int j = 1; j < M; j++)
    {
        for (int i = 0; i < N; i++) 
        {
            double y_L = A * cos(k * (i - c*j) * dx); // linke Welle
            double y_R = A * cos(k * (i + c*j) * dx); // rechte Welle
            y[i][j] = y_L + y_R; // Überlagerung
        }
    }

    // Ausgabe der Amplitude am letzten Zeitpunkt

    t = (M-1) * dt;

    for (int i = 0; i < N; i++) 
    {
        x = i * dx;
        printf("x=%f t=%f y=%f\n", x, t, y[i][M-1]);
    }

    return 0;
} 


Das Programm kann mit "gcc programm_name -lm" kompiliert werden und mit "./a.out" aufgerufen werden. Wenn man t ändert, dann man jeden beliebiegen Zeitpunkt zwischen 0 und M ausgeben. Oder man kann Amplidude Spalten mit y(t) ausgeben (Matrix).]

PS: PHP-Code sollte eigentlich c-Code heißen.

Statt die Wellen als COS-Funktionen darzustellen kann man diese auf Differenzahlgleichungen reduzieren:

PHP-Code:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define N 100 // Anzahl der Diskretisierungspunkte im Raum
#define M 1000 // Anzahl der Diskretisierungspunkte in der Zeit
#define A 1.0 // Amplitude der Welle
#define k 1.0 // Wellenzahl
#define v 1.0 // Geschwindigkeit der Welle
#define dt 0.001 // zeitliche Schrittweite
#define dx 1.0 / N // räumliche Schrittweite
#define alpha v*dt/dx // CFL-Zahl

int main()
{
    double y[N+1][M+1]; // Amplitude an den Diskretisierungspunkten
    double t, x;

    // Initialisierung der Amplitude am Anfang der Simulation

    for (int i = 0; i <= N; i++) 
   {
        y[i][0] = A * cos(k * i * dx);
    }

    // Numerische Lösung der Differentialgleichung

    for (int j = 0; j < M; j++)
    {
        // Randbedingungen setzen
        y[0][j+1] = 0.0;
        y[N][j+1] = 0.0;
        // Innere Punkte berechnen

        for (int i = 1; i < N; i++) 
       {
            y[i][j+1] = 2.0*(1-alpha*alpha)*y[i][j] - y[i][j-1] + alpha*alpha*(y[i+1][j] + y[i-1][j]);
        }
    }

    // Ausgabe der Amplitude am letzten Zeitpunkt

    t = M*dt;

    for (int i = 0; i <= N; i++) 
    {
        x = i * dx;
        printf("%f %f %f\n", x, t, y[i][M]);
    }

    return 0;
}