|
|
|
Hilfe bei Reaktions-Diffusions-Modellierung |
|
| Anna_77 |
Forum-Newbie
|
 |
Beiträge: 1
|
|
|
|
Anmeldedatum: 26.07.20
|
|
|
|
Wohnort: ---
|
|
|
|
Version: ---
|
|
|
|
 |
|
Verfasst am: 26.07.2020, 20:12
Titel: Hilfe bei Reaktions-Diffusions-Modellierung
|
 |
| |
 |
 |
Hi,
ich versuche gerade ein Diffusions-Reaktions-Gleichungssystem mit vier Größen zu programmieren. Der diskretisierte Raum ist ein Kreis mit Radius r0, der aus dem ersten System (r(1) bis r_p) und dem zweiten System (rp+1 bis r0) besteht. Im ersten System liegt der Stoff mit jeweiligen Konzentrationen vor, der sich zersetzt und die Einzelbausteine des Stoffes diffundieren durch den Stoff und schließlich in das umgebende System. Der nicht zersetzte Stoff kann nicht aus dem System entweichen.
Nun sollten diese Einzelbausteine allerdings nicht verschwinden, sondern sich in diesem äußeren System ansammeln. So, wie der Code aktuell ist, geht die Konzentration allerdings wieder auf Null zurück.
Meine Matrix ist folgendermaßen aufgebaut:
y(1) bis y(rp), y(rp +1) bis y(nr) ist die Konzentration der Bindungen in dem Stoff respektive in dem umgebenden System
y(nr+1) bis y(nr+rp), y(nr+rp +1) bis y(nr+nr) ist die Konzentration der endständigen Bindungen in dem ersten bzw. zweiten System
y(2nr+1) bis y(2nr+rp), y(2nr+rp +1) bis y(2nr+nr) ist die Konzentration der Einzelbausteine im Stoff bzw. im umgebenden System
y(3nr+1) bis y(3nr+rp), y3(nr+rp +1) bis y(3nr+nr) ist die Konzentration der durch die Einzelbausteine erzeugten H+ Ionen.
[code/]
clear all;
clc
% constants for the model
n_i = 2751;
Eps_0 = 18041;
C_c = 1;
t_c = 116; %days
n = 40;
rp = n/10;
tspan = [0 1];
% Starting conditions
for i = 1:rp
y0(i) = (n_i-3)*Eps_0/n_i;
end
for i = rp+1:n
y0(i) = 0;
end
for i = n+1:n+rp
y0(i) = 2*Eps_0/n_i;
end
for i = n+rp+1:n+n
y0(i) = 0;
end
for i = 2*n+1:2*n+rp
y0(i) = 0.0;
end
for i = 2*n+rp+1:2*n+n
y0(i) = 0;
end
for i = 3*n+1:3*n+rp
y0(i) = 10^(-7.4)*1000/C_c;
end
for i = 3*n+rp+1:3*n+n
y0(i) = 10^(-7.4)*1000/C_c;
end
[t,y] = ode23s(@(t,y) degrad_mol_rad(t,y), tspan, y0);
y1 = y(:,2*n+1);
y2 = y(:,2*n+(n-1));
figure('Name','Konzentration Einzelbausteine');
plot(t,y1,t,y2);
legend('im Stoff','im umgebenden System');
function dy = degrad_mol_rad(t,y)
nr = 40;
n = 40;
rp = nr/10;
r0 = 1;
dr = r0/(nr -1);
for i = 1:nr
r(i) = (i-1)*dr;
end
drs = dr^2;
D_p = 1;
D_w = 1;
Eps_0 = 180410;
K_COOH = 1.3490e-04;
K_W = 1.008e-14;
ke = 1;
km = 1;
kmc = 1;
kec = 1;
n_i = 2751;
C_c = 1;
t_c = 1;
for i = 1:nr
if (i ==1)
dy(i)= 0.0;
dy(i+n) = 0.0;
dy(i+2*n) = 0.0;
dy(i+3*n) = 0.0;
elseif (i == n)
dy(i) = 0;
dy(i+n) = 0;
% dy(i+2*n) = 4*D_w*(y(i+2*n-1)-2*y(i+2*n))/drs;
dy(i+2*n) = 0;
dy(i+3*n) = 0.0;
elseif (1<i) && (i<=rp) % system 1
dy(i) = -3*(km+kmc*y(i+3*n))*y(i) - (ke + kec*y(i+3*n))*y(i+n)*n_i*y(i)/((n_i-3)*Eps_0/C_c);
dy(i+n) = 2*(km+kmc*y(i+3*n))*y(i) - (ke + kec*y(i+3*n))*y(i+n)*(1-(n_i*y(i))/((n_i-3)*Eps_0/C_c));
dy(i+2*n) = (ke+kec*y(i+3*n))*y(i+n) + 0.5*(ke + kec*y(i+3*n))*y(i+n)*(1-(n_i*y(i))/((n_i-3)*Eps_0/C_c))+D_p*((y(i+2*n+1)-2*y(i+2*n)+y(i+2*n-1))/drs+ 1/r(i) *(y(i+2*n+1)-y(i+2*n))/dr) + 0.5*(ke+kec*y(i+3*n))*y(i+n)*(1-(n*y(i))/((n-3)*Eps_0/C_c));
dy(i+3*n) = K_COOH/C_c*10^3*y(i+3*n)*dy(i+2*n)/(3*(y(i+3*n))^2+...
2*K_COOH/C_c*10^3*y(i+3*n)-K_W/(C_c^2)*10^6-K_COOH/C_c*10^3*y(i+2*n));
elseif i==rp+1
dy(i) = 0;
dy(i+n) = 0;
dy(i+2*n) =dy(i+2*n-1);
dy(i+3*n) = dy(i+3*n-1);
else % umgebendes system
dy(i) = 0;
dy(i+n) = 0;
dy(i+2*n) = D_w*((y(i+2*n+1)-2*y(i+2*n)+y(i+2*n-1))/drs+1/r(i) *(y(i+2*n+1)-y(i+2*n))/dr);
dy(i+3*n) = K_COOH/C_c*10^3*y(i+3*n)*dy(i+2*n)/(3*(y(i+3*n))^2+...
2*K_COOH/C_c*10^3*y(i+3*n)-K_W/(C_c^2)*10^6-K_COOH/C_c*10^3*y(i+2*n));
end
end
dy = dy';
end
[/code]
Vielleicht hat einer von Euch eine Idee dazu oder kann mir sagen, was ich übersehe und falsch mache 
|
|
|
|
|
|
|
|
|
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.
|
|
