Differentialgleichungssystem < Maple < Mathe-Software < Mathe < Vorhilfe
|
| Status: |
(Frage) beantwortet  | | Datum: | 00:29 Mi 08.10.2008 | | Autor: | ElBarto |
| Aufgabe | Für jede der beteiligten sechs Reagenzien A, B, D, E, X, Y ergibt sich eine Differentialgleichung
A'(t) = -k1*A(t)
B'(t) = -k2*B(t) * X(t)
D'(t) = -k2*B(t)*X(t)
E'(t) = k3*X(t)
X'(t) = k1*A(t)-k2*B(t)*X(t)-k3*X(t)+k4*X(t)²*Y(t)
Y'(t) = k2*B(t)*X(t)-k4*X(t)²*Y(t)
Dabei sind ki (i=1,2,3,4) chemische Reaktionskonstanten. Dieses System wird Brüsselator genannt.
Lösen Sie mit Hilfe des Euler-Verfahrens das Differentialgleichungssystem und stellen Sie die Zeitabhängigkeiten der Reagenzien E, X, Y sowie das X-Y-Phasendiagramm dar.
Parameter: k1 = 0.0015, k3 = 0.00049, k3 = 0.877, k4 = 10
Anfangsbedingungen: A(0) = 60, B(0) = 5000, D(0) = 0, E(0) = 0, X(0) = 0, Y(0) = 0
Schrittweite: [mm] \Delta [/mm] t = 0.001
Zeitintervall: [mm] 0\le [/mm] t [mm] \le [/mm] 640 |
Hallo,
ich komme mit den vielen Gleichungen nicht klar. Eine Prozedur zum Lösen eines Differentialgleichungssystem mit zwei Gleichungen habe ich schon zu Stande bekommen:
Euler:=proc(f1,f2,x0,y10,y20,xe,h,auswahl,ausgabe)
local k0,g0,x,y1,y2,i,p:
x:=x0: y1:=y10: y2:=y20: i:=0:
if ausgabe=1 then printf("x=%f y1=%f y2=%f [mm] \n",x,y1,y2): [/mm] end if:
if auswahl=1 then p[i]:=[x,y1]: end if:
if auswahl=2 then p[i]:=[x,y2]: end if:
while (xe-x)>h do
k0:=h*f1(x,y1,y2):
g0:=h*f2(x,y1,y2):
x:=x+h: i:=i+1:
y1:=y1+k0:
y2:=y2+g0:
if ausgabe=1 then printf("x=%f y1=%f y2=%f [mm] \n",x,y1,y2): [/mm] end if:
if auswahl=1 then p[i]:=[x,y1]: end if:
if auswahl=2 then p[i]:=[x,y2]: end if:
end do:
return(convert(p,list)):
end proc:
Vielleicht hat wer von euch eine Idee was ich ändern muss um mehrere Gleichungen zu lösen. Einfach analog alles auf sechs Gleichungen anpassen? Da kam bei mir dann nur Murks raus, außerdem sieht der Quellcode dann reichlich merkwürdig aus. Da muss es doch ne elegante Lösung für geben oder?
Na ja ich hoffe sehr, dass einer von euch mir weiterhelfen kann. Ich hab schon leichte Kopfschmerzen bekommen vom vielen Brainstorming und werde mich jetzt deswegen ins Bett zurückziehen.
In diesem Sinne: Gute Nacht!
Gruß Simon
P.S. Bei Verständnisschwierigkeiten oder so, bitte einfach nochmal nachfragen.
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig  | | Datum: | 09:17 Mi 08.10.2008 | | Autor: | zahllos |
Hallo,
schreibe eine Funktion die dir zum aktuellen Zustandsvektor A, B,...,Y
die Ableitung A',B',...,Y' liefert. Im Prinzip sieht der Code dieser
Funktion so aus wie deine sechs Differentialgleichungen. Dann kannst du das Eulerverfahren für ein DGL-System ein so darstellen: Zu aktuellen Zustand (den nenne ich mal Z, ein Vektor mit sechs Komponenten) berechnest du die Ableitung Z' und integrierst über die Zeit:
Z(t+h) = Z(t) + h Z'(t)
Auf diese Weise kannst du leicht andere DGL-Systeme lösen (du musst nur die Funktion anpassen), oder ein anderes Integrationsverfahren verwenden ( dann steht an Stelle von h Z' ein komplizierterer Ausdruck).
|
|
|
| |
|
hallo ElBarto,
Das Diffgleichungssystem ist eigentlich schon sehr schön
vorbereitet, um es mit dem Eulerverfahren zu behandeln.
Ich möchte hier nur skizzenhaft das Programm bzw. die
Prozedur beschreiben:
Konstanten definieren:
k1:=0.0015 ...... k4:=10
tmin:=0 : tmax:= 640 : dt:=0.001
Anfangsbedingungen festlegen:
t:=tmin : A:=60 ...... Y:=0
While t<tmax do
Zeichnen: {Prozeduraufruf: die Prozedur "Zeichnen" soll die
zu zeichnenden Graphen um je ein kleines Stück
ergänzen}
dA:=-k1*A*dt
dB:=-k2*B*X*dt
.....
.....
[mm] dY:=(k2*B*X-k4*X^2*Y)*dt
[/mm]
t:=t+dt
A:=A+dA
B:=B+dB
.....
.....
Y:=Y+dY
end do
Falls du die berechneten Daten nicht nur grafisch ausgeben,
sondern auch in Tabellenform festhalten willst, kannst du
dies z.B. in einem Aufwasch mit der zu erstellenden Prozedur
"Zeichnen" tun.
LG
Al Chwarizmi
|
|
|
| |
|
Ich habe versucht die Idee zu realisieren. Beiliegend das
Produkt meiner Bemühungen (in Top-Pascal). Es scheint
zu funktionieren, es zeigen sich Oszillationen.
Trotzdem noch eine Frage zu den Gleichungen:
Die Grösse D scheint identisch mit B und damit überflüssig
zu sein. Könntest du die ganzen Gleichungen und
Konstanten nochmals genau überprüfen ?
Nebenbei würde ich noch gerne erfahren, welche chemischen
Reaktionen (mit welchen Reagenzien) hier genau modelliert werden.
Gruß
Al Chwarizmi
![[a]](/images/paperclip.gif "Dateianhänge") Datei-Anhang
Dateianhänge:
Anhang Nr. 1 (Typ: doc) [nicht öffentlich]
|
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig | | Datum: | 01:27 Do 09.10.2008 | | Autor: | ElBarto |
Vielen Dank erstmal für die Reaktionen. Kann mich leider erst nächste Woche wieder intensiv damit beschäftigen, dann werde ich auch auf den näheren chemischen Zusammenhang eingehen.
Gruß Simon
|
|
|
| |
|
| Status: |
(Frage) beantwortet | | Datum: | 12:45 Fr 17.10.2008 | | Autor: | ElBarto |
| Aufgabe | | siehe Ausgangspost |
Hallo,
erstmal noch einmal vielen Dank für eure Antworten. Ich hab es jetzt dann doch irgendwie geschafft reproduzierbare Ergebnisse zu bekommen. Allerdings bekomme ich nur für Werte bis x=1 Ergebnisse, danach steht bei den jeweiligen y-Werte immer "Float(undefined)".
Hier ist mein Quelltext:
> f1:=(x,y1,y2,y3,y4,y5,y6)->k1*y1; f2:=(x,y1,y2,y3,y4,y5,y6)->k2*y2*y5; f3:=(x,y1,y2,y3,y4,y5,y6)->k2*y2*y5; f4:=(x,y1,y2,y3,y4,y5,y6)->-k3*y5; [mm] f5:=(x,y1,y2,y3,y4,y5,y6)->-k1*y1+k2*y2*y5+k3*y5-k4*y5^2*y6; f6:=(x,y1,y2,y3,y4,y5,y6)->-k2*y2*y5+k4*y5^2*y6;
[/mm]
> Euler:=proc(f1,f2,f3,f4,f5,f6,x00,y10,y20,y30,y40,y50,y60,xe,h,auswahl,ausgabe)
local a0,b0,c0,e0,x0,y0,x,y1,y2,y3,y4,y5,y6,i,p:
x:=x00: y1:=y10: y2:=y20: y3:=y30: y4:=y40: y5:=y50: y6:=y60: i:=0:
if ausgabe=1 then printf("x=%f y1=%f y2=%f y3=%f y4=%f y5=%f y6=%f [mm] \n",x,y1,y2,y3,y4,y5,y6): [/mm] end if:
if auswahl=1 then p[i]:=[x,y1]: end if:
if auswahl=2 then p[i]:=[x,y2]: end if:
if auswahl=3 then p[i]:=[x,y3]: end if:
if auswahl=4 then p[i]:=[x,y4]: end if:
if auswahl=5 then p[i]:=[x,y5]: end if:
if auswahl=6 then p[i]:=[x,y6]: end if:
while (xe-x)>h do
a0:=h*f1(x,y1,y2,y3,y4,y5,y6):
b0:=h*f2(x,y1,y2,y3,y4,y5,y6):
c0:=h*f3(x,y1,y2,y3,y4,y5,y6):
e0:=h*f4(x,y1,y2,y3,y4,y5,y6):
x0:=h*f5(x,y1,y2,y3,y4,y5,y6):
y0:=h*f6(x,y1,y2,y3,y4,y5,y6):
x:=x+h: i:=i+1:
y1:=y1+a0:
y2:=y2+b0:
y3:=y3+c0:
y4:=y4+e0:
y5:=y5+x0:
y6:=y6+y0:
if ausgabe=1 then printf("x=%f y1=%f y2=%f y3=%f y4=%f y5=%f y6=%f [mm] \n",x,y1,y2,y3,y4,y5,y6): [/mm] end if:
if auswahl=1 then p[i]:=[x,y1]: end if:
if auswahl=2 then p[i]:=[x,y2]: end if:
if auswahl=3 then p[i]:=[x,y3]: end if:
if auswahl=4 then p[i]:=[x,y4]: end if:
if auswahl=5 then p[i]:=[x,y5]: end if:
if auswahl=6 then p[i]:=[x,y6]: end if:
end do:
return(convert(p,list)):
end proc:
Vielleicht sieht ja wer von euch wo der Fehler sein könnte.
Finde es auch seltsam, dass ich mir von Maple über eine halbe Million Werte ausgeben lassen soll, das dauert doch ewig bis der alle berechnet hat (wegen einem Zeitbereich von 0 [mm] \le [/mm] t [mm] \le [/mm] 640 bei einer Schrittweite von [mm] \Delta [/mm] t=0,001).
Wegen der Aufgaben habe ich nochmal alles überprüft, das stimmt so wie ich es geschrieben hab. Den chemischen Zusammenhang, sofern er aus der Aufgabenstellung hervorgeht, füge ich als Bild mit an - wenn ich wüsste wo das geht?! Man kann es aber auch unter: http://de.wikipedia.org/wiki/Br%C3%BCsselator nachlesen (dort ist der Unterschied nur, dass C dem D und D dem E aus dieser Aufgabenstellung entsprechen).
Also ich hoffe, dass mir irgendwer weiterhelfen kann bei meinem Problem und bin über jede Antwort dankbar.
MfG Simon
Edit: Ich habe gerade bemerkt, dass wenn man x später starten lässt, er auch die Werte anzeigt. Also lässt man x von 1 bis 2 gehen zeigt er dort alle Werte an. Hat dann wohl eher was mit der Menge der Werte zu tun?! Sehr merkwürdig...
|
|
|
| |
|
Als ich bei meinem Programm einmal den Zeitschritt [mm] \Delta{t}
[/mm]
deutlich vergrösserte, ergaben sich auch numerische Probleme,
es gab an einer Stelle plötzlich overflow-Meldungen.
Der kleine Zeitschritt ist also für die Berechnung schon wichtig,
wenigstens in jenen Augenblicken, wo die chemischen
Reaktionen umschlagen. Genau in dieser Phase wird auch
die numerische Lösung eventuell instabil.
Du kannst allenfalls versuchen, [mm] \Delta{t} [/mm] sogar noch kleiner
zu machen, z.B. [mm] \Delta{t}=0.0003.
[/mm]
Dagegen ist es nicht unbedingt nötig, auch alle Daten wirklich
auszugeben. Es genügt vielleicht, zehn Datenpunkte pro
Sekunde Reaktionszeit auszugeben.
LG
Übrigens gäbe es für solche Differentialgleichungen, die nur
in gewissen kleinen Gebieten starke Schwankungen beschreiben,
angepasste "adaptive" numerische Lösungsverfahren.
|
|
|
|
