|
|
| Klaatu |
Forum-Newbie
|
 |
Beiträge: 7
|
|
|
|
Anmeldedatum: 22.09.09
|
|
|
|
Wohnort: ---
|
|
|
|
Version: ---
|
|
|
|
 |
|
Verfasst am: 24.09.2009, 19:19
Titel: rekursive DGL lösen?
|
 |
Hallo zusammen.
Ich bin noch immer nicht besonders bewandert was Matlab angeht daher wende ich mich mit viel Hoffnung an euch. Zu meinem Problem:
Ich möchte ein System von DGLs der Form y'=f(y', y, x), x'=f(y',y,x) lösen. Die erste Funktion ist also rekursiv und lässt sich auch nicht umstellen. Für diskrete Punkte kann ich y' mit fzero berechnen. Gibt es eine Möglichkeit die iterative Berechnung von y' in ode45 einzubauen? Meine Versuche haben da nicht viel ergeben. Oder gibt es eine andere Möglichkeit ein solches Gleichungssystem zu berechnen?
Mit freundlichen Grüßen,
Klaatu
|
|
|
|
|
|
| Harald |
Forum-Meister
|
|
Beiträge: 24.449
|
|
|
|
Anmeldedatum: 26.03.09
|
|
|
|
Wohnort: Nähe München
|
|
|
|
Version: ab 2017b
|
|
|
|
|
|
|
Verfasst am: 24.09.2009, 21:22
Titel:
|
|
Hallo,
es sollte eigtl kein Problem sein... An ode45 oder ähnliches eine Funktion übergeben, die in etwa so aussieht:
Grüße,
Harald
|
|
|
|
| Klaatu |
Themenstarter
Forum-Newbie
|
|
Beiträge: 7
|
|
|
|
Anmeldedatum: 22.09.09
|
|
|
|
Wohnort: ---
|
|
|
|
Version: ---
|
|
|
|
|
|
|
Verfasst am: 24.09.2009, 23:07
Titel:
|
|
Hallo Harald,
vielen Dank für die erneute Hilfe. Habe das mal versucht mit der Funktion fzero. Allerdings habe ich noch Schwierigkeiten die z(1).....z(4) der Funktion zu übergeben. Des weiteren hatte ich bei kleinen Proberechnungen immer den Fehler des überschrittenen Recursionslimits. Auch eine Erhöhung hatte da nichts bewirkt außer, dass bei 5000 Matlab abschmiert. Hast du vll noch einen Tipp für mich?
Mit freundlichen Grüßen,
Klaatu
|
|
|
|
| Nicolas S. |
Forum-Century
|
|
Beiträge: 143
|
|
|
|
Anmeldedatum: 15.07.09
|
|
|
|
Wohnort: ---
|
|
|
|
Version: R2014a/b
|
|
|
|
|
|
|
Verfasst am: 25.09.2009, 08:30
Titel:
|
|
rekursiv != implizit
Ansonsten: Notfalls mal nachschauen, woran das System scheitert. Einfach mal einen impliziten Euler von Hand unter Matlab programmieren und schauen, wo es schiefgeht.
Grüße
Nicolas
_________________
--
The programmer suggested it.
|
|
|
|
| Harald |
Forum-Meister
|
|
Beiträge: 24.449
|
|
|
|
Anmeldedatum: 26.03.09
|
|
|
|
Wohnort: Nähe München
|
|
|
|
Version: ab 2017b
|
|
|
|
|
|
|
Verfasst am: 25.09.2009, 09:08
Titel:
|
|
Hallo Klaatu,
poste doch mal deinen Code bzw. das genaue Problem.
Grüße,
Harald
|
|
|
|
| Klaatu |
Themenstarter
Forum-Newbie
|
|
Beiträge: 7
|
|
|
|
Anmeldedatum: 22.09.09
|
|
|
|
Wohnort: ---
|
|
|
|
Version: ---
|
|
|
|
|
|
|
Verfasst am: 25.09.2009, 13:07
Titel:
|
|
| |
 |
 |
Hallo zusammen,
erstmal danke dafür, dass ihr euch mit meinem Problem beschäftigt.
Ich habe nochmal die Sache mit ode45 in Kombination mit fzero bearbeitet. Für einfache Systeme funktioniert das auch ganz prima. (Ich hatte den Code nochmal neu geschrieben, danach klappte es.) Ich hatte zu Testzwecken die Nullstelle einer linearen Funktion suche lassen. Eine quadratische Funktion funktionierte schon nicht mehr. Ich schätze weil fzero zwei Nullstellen zurückgibt. Für mein eigentliches Problem funktioniert es ebenfalls nicht. Die Nullstellensuche gibt zuverlässig einen Wert zurück, in Kombination mit ode45 ergeben sich allerdings konstante Werte als Lösung. Ich hab euch mal mein DGL-System als PDF angehängt. Die Terme sind recht lang. Mein Code ist recht einfach (Bis auf die Formeln) und folgender:
[Code]%DGL-funktion
function dy = KKerweitertDGLtest( t,y )
global Lp sigma omega R T d delta1 delta2 D1 D2
delta1=0.01;
delta2=0.01;
D1=5e-10;
D2=5e-10;
R=8.315;%allgemeine Gaskonstante [J/molK]
Lp=5e-12;%Parameter Lp [m³/Ns]
sigma=0.5;%Parameter sigma [-]
omega=2e-10;%Parameter omega [mol/Ns]
T=298.15;%Temperatur des Versuchs [K]
d=0.005;%Durchmesser der Membran [m]
A=0.25*pi*d^2;%Fläche der Membran
%y(1)=vfeed
%y(2)=vfermenter
%y(3)=nfeed
%y(4)=nfermenter
%c2=cfeed=y(3)/y(1)
%c1=cfermenter=y(4)/y(2)
dy = zeros(4,1); % Zeilenvektor
f=@(Jv) Lp*sigma*R*T*(exp((1-sigma)*Jv/(omega*R*T))-1)*((y(4)/y(2)-y(3)/y(1))*exp(Jv*delta2/(D2))*(1-sigma)+sigma*(y(4)/y(2)*exp(Jv*delta1/(D1))*exp(Jv*delta2/(D2))-y(3)/y(1)))/((1-sigma)*(1-exp((1-sigma)*Jv/(omega*R*T))*exp(Jv*delta1/(D1))*exp(Jv*delta2/(D2)))+exp(Jv*delta1/(D1))*sigma*(1-exp((1-sigma)*Jv/(omega*R*T))))-Jv;
z=fzero(f, 1e-12);
dy(1)=-A*z;
dy(2)=A*z;
dy(3)=-A*z*(1-sigma)*(y(4)/y(2)-y(3)/y(1)*exp((1-sigma)*z/(omega*R*T))*exp(z*delta1/D1)*exp(z*delta2/D2))/((1-sigma)*(1-exp((1-sigma)*z/(omega*R*T))*exp(z*delta1/D1)*exp(z*delta2/D2))+exp(z*delta1/D1)*sigma*(1-exp((1-sigma)*z/(omega*R*T))));
dy(4)=A*z*(1-sigma)*(y(4)/y(2)-y(3)/y(1)*exp((1-sigma)*z/(omega*R*T))*exp(z*delta1/D1)*exp(z*delta2/D2))/((1-sigma)*(1-exp((1-sigma)*z/(omega*R*T))*exp(z*delta1/D1)*exp(z*delta2/D2))+exp(z*delta1/D1)*sigma*(1-exp((1-sigma)*z/(omega*R*T))));
end
[\code]
Diese Funktion übergebe ich ode45. Als Lösung erhalte ich Konstante Werte von y(1),y(2) und y(3). Das ist nicht das was man erwarten würde, daher denke ich steckt hier irgendwo ein Fehler. Vll hat ja irgendjemand eine Idee wie ich das zum laufen bekomme.
Mit freundlichen Grüßen,
Klaatu
| Beschreibung: |
|
Download |
| Dateiname: |
DGL-System.pdf |
| Dateigröße: |
45.18 KB |
| Heruntergeladen: |
578 mal |
|
|
|
|
| Harald |
Forum-Meister
|
|
Beiträge: 24.449
|
|
|
|
Anmeldedatum: 26.03.09
|
|
|
|
Wohnort: Nähe München
|
|
|
|
Version: ab 2017b
|
|
|
|
|
|
|
Verfasst am: 25.09.2009, 14:28
Titel:
|
|
Das ist auf Anhieb schwierig zu beantworten.
Hinweise:
- die Reihenfolge der Gleichungen stimmt nicht mit dem PDF überein (Vorzeichen)
- Definiere doch A0, A1, A2 als function handles, dann wird es übersichtlicher.
- Teste mal, ob deine Funktion für t = 0 und den Anfangswert, den du übergibst, eine sinnvolle (das musst du als Anwender wissen) Rückgabe liefert
- wenn die Lösung konstant zu sein scheint, hast du sie vielleicht über einen zu geringen Zeitraum simuliert?
Grüße,
Harald
|
|
|
|
| Klaatu |
Themenstarter
Forum-Newbie
|
|
Beiträge: 7
|
|
|
|
Anmeldedatum: 22.09.09
|
|
|
|
Wohnort: ---
|
|
|
|
Version: ---
|
|
|
|
|
|
|
Verfasst am: 25.09.2009, 15:17
Titel:
|
|
Kommando zurück. Peinlicher wirds nimmer. Wenn man mit m³ rechnet und mL erwartet ist es kein wunder, dass die Werte scheinbar stagnieren. Das Volumen war um den Faktor 10^6 zu groß.
Eure Tipps waren ganz richtig und funktionieren vollkommen. Tut mir wirklich leid, dass ich hier so einen Quatsch poste. Vielen Dank für die schnelle und gute Hilfe.
Da ich noch die Versuchsdaten einlesen und die Parameter fitten muss, komme ich, wenn möglich, bestimmt noch mal darauf zurück. Diesmal dann hoffentlich mit qualifizierteren Fragen.
Gruß und besten Dank,
Klaatu
|
|
|
|
| Harald |
Forum-Meister
|
|
Beiträge: 24.449
|
|
|
|
Anmeldedatum: 26.03.09
|
|
|
|
Wohnort: Nähe München
|
|
|
|
Version: ab 2017b
|
|
|
|
|
|
|
Verfasst am: 25.09.2009, 15:22
Titel:
|
|
Ich hatte es mir ja vorhin noch verkniffen, aber in diesem Sinne dann: SI-Einheiten!!  Wenn es dich beruhigt: den meisten von uns wird sowas schon passiert sein, mir jedenfalls ja.
Grüße,
Harald
|
|
|
|
|
|
|
Einstellungen und Berechtigungen
|
|
Du kannst Beiträge in dieses Forum schreiben.
Du kannst auf Beiträge in diesem Forum antworten.
Du kannst deine Beiträge in diesem Forum nicht bearbeiten.
Du kannst deine Beiträge in diesem Forum nicht löschen.
Du kannst an Umfragen in diesem Forum nicht mitmachen.
Du kannst Dateien in diesem Forum posten
Du kannst Dateien in diesem Forum herunterladen
|
|
Impressum
| Nutzungsbedingungen
| Datenschutz
| FAQ
|  RSS
Hosted by:
Copyright © 2007 - 2024
goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks
MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.
|
|
