WICHTIG: Der Betrieb von goMatlab.de wird privat finanziert fortgesetzt. - Mehr Infos...

Mein MATLAB Forum - goMatlab.de

Mein MATLAB Forum

 
Gast > Registrieren       Autologin?   

Partner:




Forum
      Option
[Erweitert]
  • Diese Seite per Mail weiterempfehlen
     


Gehe zu:  
Neues Thema eröffnen Neue Antwort erstellen

Lösung von linearer Differentialgleichung berechnen

 

Thomas84
Forum-Meister

Forum-Meister


Beiträge: 546
Anmeldedatum: 10.02.10
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 11.02.2011, 09:54     Titel: Lösung von linearer Differentialgleichung berechnen
  Antworten mit Zitat      
Guten Tag,

ich möchte die Lösung der Differentialgleichung:

 \dot{x}(t) = f(t) \cdot x(t) + g(t) \quad \forall t

bestimmen. Die analytische Lösung dafür lautet:

 x(t) = \tilde{x}(t) (x_0 + \int_0^t \frac{g(t')}{\tilde{x}(t')} dt' ) \qquad , \tilde{x}(t) = e^{-\int_0^t f(t') dt'} \qquad , x_0 = x(0)

Für gegebene Vektoren f und g und Anfangszustand x_0 würde ich gern die Lösung x berechnen. Dabei ergeben sich einige Probleme.
 \int f dt darf nicht größer als 700 werden,da e^700 = inf. Daher muss man die Lösung gegebenenfalls zusammenstückeln. Außerdem sollte sich im Spezialfall  f= const. die gleiche Lösung wie bei lsim ergeben.

Hier ist mein bisheriger Code:

Code:

function x = lode(f,g,t,x1)
% x = lode(f,g,t,x(0))
% Lösung der linearen Differentialgleichung dx/dt = f(t)*x + g(t)
%
% x(t) = (int_0^t(g(t')/z(t') dt') + x(0))z(t)    z(t) = exp(int_0^t(f(t')dt'))
%
% Beispiel: t = linspace(0,100,1000)';
%           f = -0.1;
%           g = 100*ones(size(t));
%           x0 = 1;
%           y_lsim = lsim(f,1,1,0,g,t,x0);
%           y_lode = lode(f,g,t,x0);
%           % für konstantes f sollte sich für lode und lsim das gleiche
%           % resultat ergeben
%           figure('name','Vergleich lsim und lode');
%           subplot(2,1,1);
%           plot(t,[y_lsim,y_lode]);
%           subplot(2,1,2);
%           plot(t,y_lsim-y_lode);
%
%           figure('name','lode_bsp');
%           f = -0.1 + 0.01*sin(0.1*t);
%           y_lode = lode(f,g,t,x0);
%           plot(t,[y_lode,y_lsim]);

t = t(:);f = f(:);g = g(:);

x = zeros(size(t));
x(1) = x1;
if length(f) == 1
    f = f*ones(size(t));
end

if length(g) == 1
    g = g*ones(size(t));
end

dt = [0;diff(t)];

F = cumsum(f.*dt);
int_ss = interp1(F,1:length(t),-500:-500:F(end));
index_t = [1, floor(int_ss),length(t)];

for k = 2:length(index_t)
    i1 = index_t(k-1);
    i2 = index_t(k);
    xtilde = exp(cumsum(f(i1:i2).*dt(i1:i2)));
    korr = (exp(-f(i1:i2).*dt(i1:i2)) - 1)./(-f(i1:i2).*dt(i1:i2));
   
    integrand = g(i1:i2)./(xtilde);
    %integral = -[0;cumsum(integrand(2:end)) - cumsum(integrand(1:end-1))]./f(i1:i2);
    integral = [0;cumsum(integrand(2:end))].*dt(i1:i2);
    x(i1:i2) = (integral+ x(i1)).*xtilde./korr;
end

x(1) = x1;

 


Leider stimmen die Lösungen von lsim und lode noch nicht ganz überein (siehe Beispiel).

Hat jemand eine Idee wie man den Code verbessern könnte, oder kennt jemand eine Funktion mit der man die DGL schnell (also kein odexx) lösen kann?

viele Grüße
Thomas
Private Nachricht senden Benutzer-Profile anzeigen


Harald
Forum-Meister

Forum-Meister


Beiträge: 24.432
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 11.02.2011, 10:24     Titel:
  Antworten mit Zitat      
Hallo,

Idee wäre stichpunktartig:
- x~(t) in die Klammer ziehen
- x~(t) auch in das Integral ziehen
- Quotient zweier e-Funktionen ist e^(Differenz)
- Die dann entstehende Differenz der Integrale dürfte darauf hinauslaufen, dass du nur einmal und über ein anderes Intervall integrierst.

Grüße,
Harald
_________________

1.) Ask MATLAB Documentation
2.) Search gomatlab.de, google.de or MATLAB Answers
3.) Ask Technical Support of MathWorks
4.) Go mad, your problem is unsolvable ;)
Private Nachricht senden Benutzer-Profile anzeigen
 
Thomas84
Themenstarter

Forum-Meister

Forum-Meister


Beiträge: 546
Anmeldedatum: 10.02.10
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 12.02.2011, 07:12     Titel:
  Antworten mit Zitat      
Die Idee hatte ich auch schon. Allerdings weiss ich nicht wie ich den hinteren Integralausdruck effektiv berechnen soll.


<br />
x(t) = x_0 e^{\int_0 ^t f(t') dt'} + \int_0^t g(t') e^{\int_{t'}^t f(t'') dt''} dt'
<br />


<br />
x = test
<br />

scheinbar hat das forum oder ich heute Schwierigkeiten mit dem math-mode.
Private Nachricht senden Benutzer-Profile anzeigen
 
Harald
Forum-Meister

Forum-Meister


Beiträge: 24.432
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 12.02.2011, 10:38     Titel:
  Antworten mit Zitat      
Hallo,

zwischen der Formel und den Tags dürfen wohl keine Leerzeilen sein. Der schöneren Ansicht halber:
x(t) = x_0 e^{\int_0 ^t f(t') dt'} + \int_0^t g(t') e^{\int_{t'}^t f(t'') dt''} dt'

Zur weiteren Vorgehensweise: du hattest da doch schon einen Ansatz, und die Umformung der Gleichung sollte ihn doch eher vereinfachen?

Grüße,
Harald
_________________

1.) Ask MATLAB Documentation
2.) Search gomatlab.de, google.de or MATLAB Answers
3.) Ask Technical Support of MathWorks
4.) Go mad, your problem is unsolvable ;)
Private Nachricht senden Benutzer-Profile anzeigen
 
Thomas84
Themenstarter

Forum-Meister

Forum-Meister


Beiträge: 546
Anmeldedatum: 10.02.10
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 13.02.2011, 06:50     Titel:
  Antworten mit Zitat      
In der zweiten Form taucht im Gegensart zu meiner zuerst geposteten Lösung das Integral \int_{t'}^{t} f(t'') dt'' auf (Exponent der e-Funktion). Da dieser Ausdruck von t' und t abhängt müsste ich es als 2-dimensionalle Matrix speichern (oder eine Schleife verwenden). Die Erstellung dieser Matrix ist sicherlich sehr zeitaufwendig und sprengt auch schnell den Speicher.

Ich hatte ursprünglich eher an eine Lösung mit Hilfe von Laplace-Transformation oder ähnlichen gedacht. Da aber f zeitabhängig ist geht das wohl nicht.

viele Grüße
Thomas
Private Nachricht senden Benutzer-Profile anzeigen
 
Harald
Forum-Meister

Forum-Meister


Beiträge: 24.432
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 13.02.2011, 10:14     Titel:
  Antworten mit Zitat      
Hallo,

wenn t für den Moment mal als fest angenommen wird, sollten die verschiedenen Integrale von t' bis t als Integrale von t - k*dt bis t darstellbar sein. Ich würde es also mit cumsum/cumtrapz über den umgedrehten Vektor f(end:-1:1) versuchen.

In der nächsten Iteration der Schleife (t+dt statt t) müsstest du an den Vektor nur vorne eine 0 anfügen (für ein neues Teilstück) und zu dem Vektor als ganzen das Integral von t bis t + dt (was ja sehr schnell berechnet sein sollte) dazuzählen.

Aber mal ne ganz andere Frage: fehlt in deinem ursprünglichen Code nicht ein Minus? xtilde = exp(- integral).

Empfehlung wäre auch, den Code an sich zu dokumentieren. Ich habe da gerade meine Schwierigkeiten, das nachzuvollziehen, und fürchte, dass du selbst nach einem oder mehreren Monaten da auch Probleme haben wirst.

Grüße,
Harald
_________________

1.) Ask MATLAB Documentation
2.) Search gomatlab.de, google.de or MATLAB Answers
3.) Ask Technical Support of MathWorks
4.) Go mad, your problem is unsolvable ;)
Private Nachricht senden Benutzer-Profile anzeigen
 
Neues Thema eröffnen Neue Antwort erstellen



Einstellungen und Berechtigungen
Beiträge der letzten Zeit anzeigen:

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 | goMatlab RSS Button 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.