[ghostwhisperer]

Hi there!
Ich komm mit meiner Simulation elastischer Stöße einfach nicht weiter..
Weiß nicht wo ich das am besten reinstelle..
Im großen und ganzen stimmt mein Algo in MATLAB. Ich definiere den Flug einfacher Teilchen und ihre Stöße miteinander. Dabei definiere ich Stöße quasi so, als ob eine Feder mit der Feder-Konstante D zwischen den stoßenden Teilchen ist, sobald der Abstand kleiner einer Konstanten 2*ru ist (quasi die Radien der Kugeln).
Die Art der Abstoßung erscheint mir dabei richtig.
Wenn ich aber die Energien berechne, die alle Teilchen zu jedem Zeitpunkt haben, so kommen diese nur Betragsmäßig richtig raus.. Das Problem ist, dass kinetische und potentielle Energien, egal was ich mache, immer um einen Zeitschritt dt auseinander liegen, so dass ihre Summe nicht immer eine Konstante ist.
Frage: Was mache ich bei der Berechnung der Stöße oder den daraus folgenden Änderungen der Positionen und Geschwindigkeiten der Teilchen falsch??

Rmax = 18.0;
ru = 0.5;
D = 25;

x(1) = -1;
y(1) = 0;
z(1) = 0;
vx(1) = 1.0;
vy(1) = 0;
vz(1) = 0;
x(2) = 1;
y(2) = 0.0;
z(2) = 0;
vx(2) = -1.0;
vy(2) = 0;
vz(2) = 0;
ax(1) = 0;
ay(1) = 0;
az(1) = 0;
ax(2) = 0;
ay(2) = 0;
az(2) = 0;
dt = 0.005;
TN = 25;
Epn = 0;
ki =1;
for t = 0:dt:TN
    
    Epn = 0;
    for i = 1:AN 
        p1(1) = x(i);
        p1(2) = y(i);
        p1(3) = z(i);
        a(1) = 0;
        a(2) = 0;
        a(3) = 0;
        Fd(1) = 0;
        Fd(2) = 0;
        Fd(3) = 0;
        
        for j = 1:AN
            if i ~= j
                p2(1) = x(j);
                p2(2) = y(j);
                p2(3) = z(j);
                if norm(p2-p1,2)<(2*ru)
                    s = p2 - p1;
                    sn = s/norm(s,2);
                    s = s - 2*ru*sn;
                  if norm(s,2)<(2*ru-0.02*ru)
                        s = s;
                    else
                        s =0*s;
                    end
                    Fd = Fd -D*s;
                    Epn = Epn + 0.5*D*norm(s,2)*norm(s,2);
                end
            end    
        end

       a = -Fd / m;
        
        ax(i) = a(1);
        ay(i) = a(2);
        az(i) = a(3); 
    end
    
 Ep(ki) = 0.5*Epn;

    for i = 1:AN 
        
        p1(1) = x(i);
        p1(2) = y(i);
        p1(3) = z(i);
        v1(1) = vx(i);
        v1(2) = vy(i);
        v1(3) = vz(i);
        a(1) = ax(i);
        a(2) = ay(i);
        a(3) = az(i);
        
        if abs(p1(1))>=Rmax || abs(p1(2))>=Rmax || abs(p1(3))>=Rmax
            if abs(p1(1))>=Rmax 
                if (v1(1)/abs(v1(1))) == (p1(1)/abs(p1(1)))
                    v1(1) = -v1(1)
                end
            end
            if abs(p1(2))>=Rmax 
                if (v1(2)/abs(v1(2))) == (p1(2)/abs(p1(2)))
                    v1(2) = -v1(2)
                end
            end
            if abs(p1(3))>=Rmax 
                if (v1(3)/abs(v1(3))) == (p1(3)/abs(p1(3)))
                    v1(3) = -v1(3)
                end
            end
        end
        
        Ek(ki) = Ek(ki)+ 0.5*(v1(1)*v1(1)+v1(2)*v1(2)+v1(3)*v1(3));
        dv = a*dt;
        v1 = v1 + dv;
        p1 = p1 + v1*dt;
        
        vx(i) = v1(1);
        vy(i) = v1(2);
        vz(i) = v1(3);
        x(i) = p1(1);
        y(i) = p1(2);
        z(i) = p1(3); 
    end
      
    ti(ki) = t;
   ki=ki+1;

end
plot(ti,Ek+Ep);
plot(ti,Ek,'-',ti,Ep,'--',ti,Ek+Ep,'.');
f=vq;
end

PS: Initialisierungen und unrelevanten Code hab ich nicht hierherkopiert.

DANKE!