Lösung von linearer Differentialgleichung berechnen

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

Partner:

Forum

Option

[Erweitert]

• Diese Seite per Mail weiterempfehlen

Gehe zu:

Thomas84

Forum-Meister


	Beiträge: 546

	Anmeldedatum: 10.02.10

	Wohnort: ---

	Version: ---

Verfasst am: 11.02.2011, 09:54 Titel: Lösung von linearer Differentialgleichung berechnen

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;

Funktion ohne Link?

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

Harald

Forum-Meister


	Beiträge: 24.495

	Anmeldedatum: 26.03.09

	Wohnort: Nähe München

	Version: ab 2017b

Verfasst am: 11.02.2011, 10:24 Titel:

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 ;)

Thomas84

Themenstarter

Forum-Meister


	Beiträge: 546

	Anmeldedatum: 10.02.10

	Wohnort: ---

	Version: ---

Verfasst am: 12.02.2011, 07:12 Titel:

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.

Harald

Forum-Meister


	Beiträge: 24.495

	Anmeldedatum: 26.03.09

	Wohnort: Nähe München

	Version: ab 2017b

Verfasst am: 12.02.2011, 10:38 Titel:

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 ;)

Thomas84

Themenstarter

Forum-Meister


	Beiträge: 546

	Anmeldedatum: 10.02.10

	Wohnort: ---

	Version: ---

Verfasst am: 13.02.2011, 06:50 Titel:

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

Harald

Forum-Meister


	Beiträge: 24.495

	Anmeldedatum: 26.03.09

	Wohnort: Nähe München

	Version: ab 2017b

Verfasst am: 13.02.2011, 10:14 Titel:

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 ;)


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 |

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.