dx/dy beim lösen von ode mit ausgeben und beim lösen nutze

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

Partner:

Forum

Option

[Erweitert]

• Diese Seite per Mail weiterempfehlen

Gehe zu:

RonG

Forum-Newbie


	Beiträge: 3

	Anmeldedatum: 15.05.14

	Wohnort: Singapore

	Version: ---

Verfasst am: 15.05.2014, 08:11 Titel: dx/dy beim lösen von ode mit ausgeben und beim lösen nutze

Hallo, nachdem ich bereits Wochen an diesem Problem verzweifle und mir auch google noch nicht weitergeholfen hat Embarassed

, versuche ich mein glück mal hier.

Mein Problem ist folgendes:

Beim Lösen von mehreren DGLs mit ode45 will ich mir neben dem Wert für x auch den Wert für dx/dy mit in die selbe Matrix schreiben lassen.

x ist bei mir hierbei die Temperatur und y die Zeit bei einer ablaufenden Reaktion 1. Ordnung.

Die Wärmekapazität welche ich zum berechnen der Heizrate brauche hängt nun genau von dieser ab...

Welche Möglichkeiten gibt es dieses Problem zu lösen und sich alle Ergebnisse in eine zeitabhängige Matrix schreiben zu lassen?

Hier mal meine Skripte für das vereinfachte Program und die DGLs:

1. main program

Code:

clear
clc

%kinetic

EA = 147.3; %kJ mol-1 Aktivierungsenergie
k0 = 60*5E15; %min-1 Präexponentialfaktor
H = 1290; %J g-1 Reaktionsenthalpie
R = 8.314E-3; %kJ/mol-1/K-1 Gaskonst.

%sample parameters

m_dtbp = 326.682E-3; %g mass dtbp for experiment
cp_dtbp = 1.8; %kJ kg-1 K-1 Spezifische Wärme der Reaktionsmischung

C_dtbp = m_dtbp*cp_dtbp; %J/K

%sample container

mTi = 844E-3; %g
cpTi = 0.561; %J/(g*K)
mSteel = 22.7; %g
cpSteel= 0.466; %J/(g*K)

C_Ti = cpTi*mTi; %J/K
CSteel = cpSteel*mSteel; %J/K

Cmin = C_Ti+C_dtbp;

%operation parameters
t1 = 0; %min start time
x1 = 0.001; %- Anfangsumsatz
T1 = 378.15; %K Anfangstemperatur (105°C)
HR1 = 0; %K/min Anfangsheizrate

%_____________________________________________________________________

%Zusammengefasst

Z1=-EA/R;
Z2=H*m_dtbp/(Cmin+CSteel*0.5); %K %H*m/C

%Anfangswerte
C1=[t1 x1 T1 HR1];

%Berechnungen

[Time Ergebnisse]=ode45(@(t,C) dgl_dtbp_test(t,k0,C,Z1,Z2),[0 200],C1); %1

%_____________________________________________________________________

%Ableitung dT/dt

i = 1;

while(i <= length(Ergebnisse(:,3)))
heat1(i) = Z2*k0*exp(Z1/Ergebnisse(i,3))*(1-(Ergebnisse(i,2)));
i = i + 1;
end

%_____________________________________________________________________

Tnew=Ergebnisse(:,3)-273.15; %Temperatur in °C

%dT/dt

subplot(2,2,1),semilogy(Tnew,heat1);
xlabel('temperature T in °C')
ylabel('dT/dt in K/min')
title('heat rate','FontSize',12)

%temperature

subplot(2,2,3),plot(Time,Tnew);
xlabel('time t in min')
ylabel('temperature T in K')
title('temperature','FontSize',12)

%turnover

subplot(2,2,2),plot(Time,1-Ergebnisse(:,2));
xlabel('time t in min')
ylabel('1-x')
title('turnover','FontSize',12)

%new HR

subplot(2,2,4),semilogy(Tnew,Ergebnisse(:,4));
xlabel('temperature T in °C')
ylabel('dT/dt in K/min')
title('heat rate - not working','FontSize',12)

%Hintergrundfarbe und Fenstergröße

set(gcf, 'color', 'white');
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);

Funktion ohne Link?

2. DGLs

Code:

function x = dgl_dtbp_test(t,k0,C,Z1,Z2)

x=zeros(4,1);

x(1) = 1; %time in min
x(2) = k0*exp(Z1/C(3))*(1-(C(2))); %turnover
x(3) = Z2*x(2); %temperature
x(4) = x(3); %heat rate - not working

end

%summarized

%Z1=-EA/R;
%Z2= H*m/C;

Funktion ohne Link?

Schonmal vielen Dank im Voraus für eure Vorschläge!

Jan S

Moderator


	Beiträge: 11.057

	Anmeldedatum: 08.07.10

	Wohnort: Heidelberg

	Version: 2009a, 2016b

Verfasst am: 17.05.2014, 13:06 Titel: Re: dx/dy beim lösen von ode mit ausgeben und beim lösen n

Hallo RonG,

ODE45 gibt die Ableitung nicht direkt zurück, aber Du hast ja die Funktion, die sie ausrechnet:

Code:

[Time, Ergebnisse] = ode45(@(t,C) dgl_dtbp_test(t,k0,C,Z1,Z2),[0 200],C1);
dxdt = dgl_dtbp_test(Time, k0, Ergebnisse.', ,Z1, Z2);

Funktion ohne Link?

Die Funktion muss dann noch Vektoren als Inputs berücksichtigen:

Code:

function x = dgl_dtbp_test(t,k0,C,Z1,Z2)
x = zeros(4, numel(t));
x(1, :) = 1; %time in min
x(2, :) = k0 * exp(Z1 ./ C(3, :)) .* (1 - C(2, :)); %turnover
x(3, :) = Z2 * x(2, :); %temperature
x(4, :) = x(3, :); %heat rate - not working
end

Funktion ohne Link?

Möglicherweise ist das mit dem Transponieren noch falsch, aber das sollte kein ernstes Problem darstellen.

Gruß, Jan

RonG

Themenstarter

Forum-Newbie


	Beiträge: 3

	Anmeldedatum: 15.05.14

	Wohnort: Singapore

	Version: ---

Verfasst am: 20.05.2014, 03:15 Titel:

Hallo Jan,

vielen Dank erstmal für deine Antwort, jedoch weiß ich noch nicht so recht, was es mir bringt die Inputs als Vektoren zu schreiben, da dies bei ode45 mit x(1) ja schon automatisch gemacht wird, oder?

Jedenfalls werden dadurch genau die selben Ergebnisse ausgegeben wie vorher auch schon.

Was ich eher bräuchte wäre folgendes, jedoch weiß ich nicht genau ob und wie man dies Matlab beibringen kann.

Code:

x=zeros(4,numel(t));

x(1, :) = 1; %time in min
x(2, :) = k0*exp(Z1/C(3, :))*(1-(C(2, :))); %turnover
x(3, :) = Z2/Cmin*x(2, :); %temperature
C(4, :) = x(3, :); %heat rate (not working)

Funktion ohne Link?


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.