Wer kennt sich in Physik gut aus ? (Runge-Kutta,ODE Solver)

[Hellspinder]

Wer kann mir eine Differenzialgleichung zur Beschleunigung eines Autos lösen ?

Phase 1 Beschleunigung bis 100KM/H
Phase 2 Leerlauf weiterrollen
Phase 3 Bremsen

Das ganze im Ode Solver nach Runge-Kutta :freak:

Ich blick da nicht durch

[Pinoccio]

Das ist nix für einen Physiker, sondern für einen Numeriker oder Ingenieur.

Könntest du die Frage etwas ausführlicher formulieren und auch mal deinen Ansatz schildern? Und meinst du mit ODE SOlver das? (http://www.javaview.de/services/odeSolver/index.html)

mfg

[BlödBär]

Ich verstehe nur Bahnhof. Ich hatte in Physik andere Sachen... ;)

[pest]

das heißt Differential
und so eine Gleichung für ein Auto aufzustellen ist ja nicht schwer, total simpel: 1. Phase f''>0, 2. Phase: f''=0, 3. Phase, f'' <0 mit f'' 2. Abl. des Weges nach der Zeit
du musst sie ja noch nichtmal selbst ausrechnen :eek:

[Pinoccio]

das heißt Differential:freak: (http://www.wissen.de/wde/generator/wissen/ressorts/bildung/woerterbuecher/index,page=2074452.html) <-Link

mfg

[pest]

:usad: (http://www.googlefight.com/index.php?lang=en_GB&word1=Differenzial&word2=Differential) <- Link

[Pinoccio]

:usad: (http://www.googlefight.com/index.php?lang=en_GB&word1=Differenzial&word2=Differential) <- LinkIch bin doch auf deiner Seite ...
Nur die amtliche Rechtschreibung (http://www.duden.de/deutsche_sprache/sprachwissen/rechtschreibung/crashkurs/fremdwortschreibung/regel_06.php) sieht eben z als Regel und t als Ausnahme.

Wie wäre: dt heißt Differential, dz heißt Differenzial und so weiter ... ;D

Und nun zurück zum Thema!

mfg

[Hellspinder]

Das ist nix für einen Physiker, sondern für einen Numeriker oder Ingenieur.

Könntest du die Frage etwas ausführlicher formulieren und auch mal deinen Ansatz schildern? Und meinst du mit ODE SOlver das? (http://www.javaview.de/services/odeSolver/index.html)

mfg

Ja sowas mein ich. Das ganze halt in Java programmieren. Das scheint soweit richtig zu sein allerdings beschleunigt das Auto erst langsam und dann bis Tempo 100 immer schneller (in der Realität eher umgekehrt).

float ODE(float gradient, float Value, float rest, float timeQuant) {
float m1, m2, m3, m4, m; // Gradienten Zwischenergebnisse
float res1, res2, res3, res; // Zwischerergebnisse
m1 = gradient*Value + rest;
res1 = Value + m1*timeQuant/2;
m2 = gradient*res1 + rest;
res2 = Value + m2*timeQuant/2;
m3 = gradient*res2 + rest;
res3 = Value +m3*timeQuant;
m4 = gradient*res3 + rest;
m = (m1+2*(m2 + m3)+m4)/6;
res = Value + m*timeQuant;
return res;
}

ich blick da nicht durch. Die Aufgabe ist :

siehe Bild

[pest]

tipp: wenn variablen temporär sind, musst du sie nicht res1,res2... nennen, erhöht die übersichtlichkeit

und ja


res = Value + m*timeQuant;


ist eine Zeitintegration

poste den Rest des Codes ;) - oder teste ob dein Verfahren funktioniert wenn du statt Runge-Kutta, Euler-Cauchy verwendest

[noid]

Wenn eine konstante Beschleunigung auf eine Körper einwirkt, dann wird ohne Luftwiderstand und Reibung das Ding immer schneller. Egal in welche Richtung der Körper beschleunigt wird.

In der Praxis sieht das natürlich anders aus, da haben wir ja kein konstantes a, und und und...
(einfach mal in einer ruhigen Minute mal drüber nachdenken)

mal zur illustration: http://www.walter-fendt.de/ph14d/beschleunigung.htm

[Hellspinder]

naja genau genommen gibt es einen Reibungskoeffizienten :

Ich hab mal den Code für die einzelnen Phasen angehängt vielleicht kann jemand den Fehler finden warum die Werte nicht passen. :confused:.
Irgendwas muss mit dem Maßstab falsch sein. Ich kriegs nich hin.
** CENSORED **

EDIT : 500 Posts !!!!!!!!!!!!!!!!!!!!!!!!!!!!
http://images.mzzt.net/awesome.png

[pest]

das sieht ziemlich falsch aus, wenn ich es in 15min nicht hinbekomme, meld ich mich morgen...

du musst natürlich vieles auf dt umrechnen und die einheiten beachten...but thats the aufgabe for pinoccio

[pest]

so habs mal schnell in C geschrieben (10 min mehr gebraucht :D)

als Zeitintegration habe ich jetzt Euler-Cauchy genommen
Aussderm ist ein Vorzeichenfehler bei der Reibungskraft Bremskraft in der Aufgabenstellung.

Einheiten müssen nach Intuition stimmen :p

Ausgabe mit Zeitschritt dt=0.01s und r=1 kg/s (keine Reibung)
die Anfangsbeschleunigung ist 100/(3.6*10)=2.78 m/s²


t: 10.0 [s] v: 99.22 [km/h] s: 138.18 [m]
t: 20.0 [s] v: 98.24 [km/h] s: 412.43 [m]
t: 22.7 [s] v: 0.18 [km/h] s: 449.45 [m]



Ausgabe mit Zeitschritt dt=0.01s und r=100 kg/s
Anfangsbeschleunigung ~4.4m/s²


t: 10.0 [s] v: 100.16 [km/h] s: 162.06 [m]
t: 20.0 [s] v: 36.83 [km/h] s: 337.80 [m]
t: 21.0 [s] v: 0.13 [km/h] s: 342.65 [m]



Quellcode



double ODE(float val,float f,float h)
{
return val + f*h;
}

double Status(float t,float v,float x)
{
printf("t: %0.1f [s] v: %0.2f [km/h] s: %0.2f [m]\n",t,v*3.6,x);
}

int main() {
float t,dt,v,m,a,r,x,F;
x = 0;
dt = 0.01;
t = 0;
m = 1000;
r = 1;
v = 0;
a = 2.77;
F = 10000;

while (t < 10) {
v = ODE(v,-r/m*v+a,dt);
t = t + dt;
x = x + dt*v;
}
Status(t,v,x);
while (t < 20) {
v = ODE(v,-r/m*v,dt);
t = t + dt;
x = x + dt*v;
}
Status(t,v,x);
while (v > 0.1) {
v = ODE(v,-r/m*v-F/m,dt);
t = t + dt;
x = x + dt*v;
}
Status(t,v,x);
return 0;
}

[Hellspinder]

unfassbar es gibt wirklich Leute die machen sowas freiwillig :freak:.

Vielen Dank für deine Mühe aber ich muss das zwingend mit Runge-Kutta machen ist das soviel anders ? und in dem Programm "Processing" das ist Java Umgebung.

[pest]

unfassbar es gibt wirklich Leute die machen sowas freiwillig :freak:.


ich habe das gerade als Entspannung vom Lernen gemacht :freak:


Vielen Dank für deine Mühe aber ich muss das zwingend mit Runge-Kutta machen ist das soviel anders ?

nö, du musst nur die Funktion ODE ersetzen, aber da muss ich nochmal drüber nachdenken wie man das geschickt (und richtig) macht.
dein Fehler lag ja schon im falschen Umformen der Ausgangsgleichung

morgen ;)

[Trap]

Wenn man es ganz nach festem Schema machen möchte, formt Differentialgleichungen höherer Ordnung erst nach http://de.wikipedia.org/wiki/Gew%C3%B6hnliche_Differentialgleichung#Reduktion_von_Gleichungen_h.C3.B6herer_Or dnung_auf_Systeme_erster_Ordnung in Differentialgleichungen erster Ordnung um.
Dann braucht man aber einen Runge-Kutta-Löser, der mit mehrdimensionalen Statusvektoren zurechtkommt.

Der Code von pest ist meiner Meinung nach falsch, weil er nur eine Ordnung/Dimension des Problems mit dem vorgegebenen Solver löst, die Integration über die Geschwindigkeit macht er immer mit Euler. Außerdem ist das Euler-Cauchy falsch implementiert, richtig wäre x_n+1=x_n+dt*v_n, der Code macht aber x_n+1=x_n+dt*v_n+1.

[pest]

Der Code von pest ist meiner Meinung nach falsch, weil er nur eine Ordnung/Dimension des Problems mit dem vorgegebenen Solver löst, die Integration über die Geschwindigkeit macht er immer mit Euler.


ich mache (bis jetzt) alles über Euler


Außerdem ist das Euler-Cauchy falsch implementiert, richtig wäre x_n+1=x_n+dt*v_n, der Code macht aber x_n+1=x_n+dt*v_n+1.


mal sehen


geg: y''=f(x,y,y') mit AB

mit z1=y und z2=y'

folgt z'1=z2=y' und z'2=f(x,z1,z2)

damit ergeben sich folgende (explizite) Integrationsgleichungen

z1_(n+1)=z1_(n) + h*z2_(n)
z2_(n+1)=z2_(n) + h*f(x_(n),z1_(n),z2_(n))


damit stimmt meine Vorgehensweise wenn man das System 2. Ordnung auf ein System gewöhnlicher DGLen reduziert.
Allerdings hattest du Recht, das ich den Weg (z1) mit dem falschen Wert berechnet habe. Kam vermutlich aus dem Gedanken, das man ja das gesamte Zeitintervall in einem Schritt integrieren
könnte, wenn man das (nach dem mathematisch korrekten Weg) macht, kommt für den Weg = 0m heraus.

[Trap]

Für die Differentialgleichungen in der Aufgabe macht das wahrscheinlich keinen großen Unterschied welche Geschwindigkeit man nimmt oder mit welcher Methode man sie zum Ort integriert.

Der Geschwindigkeitsverlauf ist ja stetig, monoton und auch sonst in jeder Hinsicht unkompliziert. x'=-1/(x+0.1) oder so wäre da spannender.

[pest]

wie schon geschrieben, der Wegintegrationsschritt war nicht ganz korrekt, aber sonst stellt die Vorgehensweise ein korrektes explizites Euler-Verfahren für ein System 2. Ordnung dar.
Ich kuck mal wie man das mit Kutta machen kann.

[Trap]

Für eine ordentliche Implementierung braucht man eine andere Programmstruktur:
interface ODE { StateVector value(StateVector a, float time);}
Statevector RK4(StateVector start, float time, ODE o, float stepsize) {
StateVector k1 = o.value(start,time);
StateVector k2 = o.value(start+stepsize/2*k1,time+stepsize/2);
StateVector k3 = o.value(start+stepsize/2*k2,time+stepsize/2);
StateVector k4 = o.value(start+stepsize*k3,time+stepsize);
return start+stepsize/6*(k1+2*k2+2*k3+k4);
}

main(){
ODE problem1=new ODE(){StateVector value(StateVector a, float time)
{StateVector result= new StateVector(2);
result[0]=a[0];
result[1]=-R_aus_Aufgabenstellung*a[1]/m+A_aus_Aufgabenstellung;
return result;};

StateVector state = new StateVector(2);
for(float t=0;t<=tEND;t+=stepsize) {
state=RK4(state, t, problem1, stepsize);
}
}


Damit es keine fertige Lösung ist hab ich mal Operatoroverloading für StateVector benutzt und StateVector nicht genauer definiert ;)

[pest]

ich will da nur i,j,x,v,a usw. sehn...grausam :D - außerdem ist deins nur für Systeme 1. Ordnung, und da fehlt noch die Bremskraft
aber du legst es ja offensichtlich drauf an, ich mach dir den minimal C-Code inkl. Euler + Kutta

[Trap]

ich will da nur i,j,x,v,a usw. sehn...grausam :D
Kein Problem: http://www.unige.ch/~hairer/prog/nonstiff/dopri5.f
außerdem ist deins nur für Systeme 1. Ordnung
Das macht nix, es gibt für jedes Problem höherer Ordnung ein dazu äquivalentes in erster Ordnung.

[pest]

Das macht nix, es gibt für jedes Problem höherer Ordnung ein dazu äquivalentes in erster Ordnung.

schon klar, wie du sicher bemerkt hast, habe ich das vor einer seite gemacht
aber es nützt dem ts wenig, allerdings sehe ich mich aufgrund deiner anwesenheit dazu gezwungen meinen code auf unter 30 Zeilen zu quetschen :D

[Migrator]

Mathematiker unter sich :freak:

http://brightsblog.files.wordpress.com/2007/05/miraclel.jpg

Ich glaub außer euch beiden peilt hier niemand mehr was :D (okay Pinoccio und so..)

[pest]

schlimm genug

[pest]

soooo :D

der Quellcode entspricht nun meinen ästhetischen Ansprüchen
allerdings glaube ich, dass ihn niemand lesen kann, aber hübsch isser

seltsam ist nur, das die Wegintegration beim RK-Verfahren extrem ungenau
wird wenn der Zeitschritt groß ist. Einen Fehler im Programm schließe ich mal aus.

Ausgabe eines Programmdurchlaufes


Zeitschritt: 0.0100 (s)

Explizites Euler-Verfahren für System 2. Ordnung
r= 1 (kg/s), a=2.778 (m/s²)
t: 10.0 [s] v: 99.50 [km/h] s: 138.29 [m]
t: 20.0 [s] v: 98.51 [km/h] s: 413.31 [m]
t: 22.7 [s] v: -0.26 [km/h] s: 450.82 [m]
r= 100 (kg/s), a=4.400 (m/s²)
t: 10.0 [s] v: 100.16 [km/h] s: 161.79 [m]
t: 20.0 [s] v: 36.83 [km/h] s: 337.70 [m]
t: 21.0 [s] v: -0.23 [km/h] s: 342.65 [m]

Explizites Runge-Kutta-Verfahren 4. Ordnung für System 2. Ordnung
r= 1 (kg/s), a=2.778 (m/s²)
t: 10.0 [s] v: 99.50 [km/h] s: 138.98 [m]
t: 20.0 [s] v: 98.51 [km/h] s: 415.38 [m]
t: 22.7 [s] v: -0.26 [km/h] s: 453.08 [m]
r= 100 (kg/s), a=4.400 (m/s²)
t: 10.0 [s] v: 100.13 [km/h] s: 162.54 [m]
t: 20.0 [s] v: 36.83 [km/h] s: 339.32 [m]
t: 21.0 [s] v: -0.21 [km/h] s: 344.30 [m]


Quellcode



#define NS 4
#define Eval(v,r,a,F) (-(r)/(m)*(v)+(a)-(F)/(m))
#define SumKT ((kt[1]+2*(kt[2]+kt[3])+kt[4])/6.)
#define SC ((state[2]+coef[i]*kt[i]))

const
float coef[4]={1,0.5,0.5,1};
float dt=0.01;
float m=1000;
float F=10000;

void Stats(float *state) {
printf(" t: %0.1f [s] v: %6.2f [km/h] s: %6.2f [m]\n",state[0],state[2]*3.6,state[1]);
}

void Euler_Step(float *state,float max_t,float a,float r,float F) {
while (state[0] < max_t) {
state[0] += dt;
state[1] += dt*state[2];
state[2] += dt*Eval(state[2],r,a,F);
if (state[2] < 0.) break;
};
Stats(state);
}

void RK_Step(float *state,float max_t,float a,float r,float F)
{
int i;
float kt[NS+1]={0};
while (state[0] < max_t) {
state[0] += dt;
for (i=0;i<NS;i++) kt[i+1] = dt*SC;
state[1] += SumKT;
for (i=0;i<NS;i++) kt[i+1] = dt*Eval(SC,r,a,F);
state[2] += SumKT;
if (state[2] < 0.) break;
};
Stats(state);
}

void Euler(float r,float a) {
float state[3]={0};
printf(" r=%4.0f (kg/s), a=%0.3f (m/s²)\n",r,a);
Euler_Step(state,10,a,r,0);
Euler_Step(state,20,0,r,0);
Euler_Step(state,30,0,r,F);
}

void RungeKutta(float r,float a) {
float state[3]={0};
printf(" r=%4.0f (kg/s), a=%0.3f (m/s²)\n",r,a);
RK_Step(state,10,a,r,0);
RK_Step(state,20,0,r,0);
RK_Step(state,30,0,r,F);
}

int main(void) {
printf("Zeitschritt: %0.4f (s)\n",dt);
printf("\nExplizites Euler-Verfahren für System 2. Ordnung\n");
Euler(1,100/(10*3.6));
Euler(100,4.4);
printf("\nExplizites Runge-Kutta-Verfahren 4. Ordnung für System 2. Ordnung\n");
RungeKutta(1,100/(10*3.6));
RungeKutta(100,4.4);
return 0;
}