Simulation EM-Feld (Matlab)

[M. M. (blackcow)]

Ich wollte mit matlab eine 1-Dimensionale EM-Welle Simulieren. Mein 
Ausgangspunkt sind die Maxwellschen Gleichungen (zur Vereinfachung habe 
ich E=D und B=H gesetzt:
rot(E)=dB/dt
rot(B)=dE/dt

Bei einer Dimension habe ich nun in Matlab sechs 1-Dimensionale Arrays 
erstellt:
buf_rotE = aktuelle Rotation von E-Feld
buf_rotB = aktuelle Rotation von B-Feld
buf_E =    vergangene&aktuelle E-Feld-Vektoren
buf_B =    vergangene&aktuelle B-Feld-Vektoren

Ein Simulationsschritt sieht bei mir folgendermaßen aus:
Aus buf_rotE und buf_rotB werden jeweils das aktuelle E- und B-Feld 
berechnet. Mit diesen Werten und den des Schrittes davor werden dann die 
Differenzen (sprich zeitliche Ableitung) berechnet und als neue 
Rotationswerte in buf_rotE, bzw. buf_rotB eingetragen.

Das Problem ist, es funktioniert nicht. Das Ergebnis sieht überhaupt 
nicht nach Welle aus. Kann mir jemand sagen wo der Haken liegt?

% Initialisiere Variablen
    % Einstellungen
        % Feldgröße
        X_width = 200;
    % Simulationsgeschwindigkeit
    d_t  = 0.01;
    % Verknüpfung rotE <- dB/dt
    k_rotE_dB = 1;
    % Verknüpfung rotB <- dE/dt
    k_rotB_dE = 1;

  buf_rotE = zeros(X_width,1);
  buf_rotB = zeros(X_width,1);
  buf_E = zeros(X_width,2);    % X, alt/neu
  buf_B = zeros(X_width,2);       % X, alt/neu
  i_alt = 1;
  i_neu = 2;
  gw_rotE  = zeros(2, 1);          % min/max
  gw_offset = 0;
  gw_skng = 0;
  X = 0;
    rain = 0;
    
buf_rotE(X_width/2) = 250;
% Führe Simulationsschritt durch
for step = 1:50
% buf_rotE(X_width/2) = 250*sin(step*0.4);
    figure(1);
    subplot(4,1,1);
        bar(buf_rotE);
    title(['step: ' num2str(step)]);
        axis([1 X_width -250 250]);
    subplot(4,1,3);
        bar(buf_rotB);
        axis([1 X_width -250 250]);
        title('rotB');
    pause(d_t);
  % Rotiere Array Index: alt <=> neu
        if i_alt==1
            i_alt = 2;
      i_neu = 1;
        else
            i_alt = 1;
            i_neu = 2;
        end
  % Berechne Felder und schreibe min/max Werte mit
        gw_rotE(1) = 0;
        gw_rotE(2) = 0;
        gw_B(1) = 0;
        gw_B(2) = 0;
        for X = 1:X_width
            % aktuelles E-Feld
                buf_E(X, i_neu) = 0;
                if X>=X_width
                    buf_E(X, i_neu) = 0; %buf_rotE(X)-buf_rotE(X-1);
                elseif X==1
                    buf_E(X, i_neu) = 0; %buf_rotE(X+1)-buf_rotE(X);
                else
                    buf_E(X, i_neu) = buf_E(X, i_alt) + buf_rotE(X)-buf_rotE(X-1);
                    buf_E(X, i_neu) = buf_rotE(X+1)-buf_rotE(X) + buf_E(X, i_neu);
                end
                % aktuelles B-Feld
                    buf_B(X, i_neu) = 0;
                    if X>=X_width
                        buf_B(X, i_neu) = 0; %buf_rotB(X)-buf_rotB(X-1);
                    elseif X==1
                        buf_B(X, i_neu) = 0; %buf_rotB(X+1)-buf_rotB(X);
                    else
                        buf_B(X, i_neu) = buf_B(X, i_neu) + buf_rotB(X+1)-buf_rotB(X);
                        buf_B(X, i_neu) = buf_rotB(X)-buf_rotB(X-1) + buf_B(X, i_neu);
                    end
        end
        for X = 1:X_width
                % aktuelles E-Rotation
                    buf_rotE(X) = ( buf_B(X, i_neu)-buf_B(X, i_alt) ) * k_rotE_dB;
                    % schreibe min/max Werte mit
                    if buf_rotE(X)<gw_rotE(1)
                        gw_rotE(1) = buf_rotE(X);
                    end
                    if buf_rotE(X)>gw_rotE(2)
                        gw_rotE(2) = buf_rotE(X);
                    end
                % aktuelles B-Rotation
                    buf_rotB(X) = ( buf_E(X, i_neu)-buf_E(X, i_alt) ) * k_rotB_dE;
        end

    subplot(4,1,2);
        bar(buf_E(:,i_neu));
        axis([1 X_width -250 250]);
        title('E');
    subplot(4,1,4);
        bar(buf_B(:,i_neu));
        axis([1 X_width -250 250]);
        title('B');
    hold off;
    % pause;
end

[Rolf R. (Gast)]

Ich verstehe nichts von Matlab.

Aber warum nimmst du nicht die Wellengleichung, die man aus den 
Maxwellschen Gleichungen erhält, und löst diese zuerst mathematisch auf 
dem Papier für deinen Anwendungsfall. Mit dem Ergebnis hast du dann dein 
E- und dein B-Feld in Abhängigkeit von Zeit und Ort.

Damit würde ich dann in Matlab weiter machen.

Das Einsetzten der Maxwellschen Gleichungen ineinander und die Umformung 
zur Wellengleichung und die Lösung dieser ist bei mir aber schon 
Jahrzehnte her. Ich müsste mich da auch wieder einarbeiten.

[Rolf R. (Gast)]

Oder willst du FDTD wie hier

https://www.youtube.com/watch?v=RoUrqR71UWo

machen?

[Bolle (Gast)]

Matlab ist dafür nicht so geeignet. Ich empfehle Comsol.