Implizites Euler-Verfahren

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

Partner:

Forum

Option

[Erweitert]

• Diese Seite per Mail weiterempfehlen

Gehe zu:

mieze579

Forum-Newbie


	Beiträge: 6

	Anmeldedatum: 22.01.15

	Wohnort: ---

	Version: ---

Verfasst am: 19.06.2015, 20:53 Titel: Implizites Euler-Verfahren

Hallo zusammen,

ich habe einige Gleichungen, die ineinandergreifen und soll diese Gleichungen mit dem IMPLIZITEN EULER-VERFAHREN lösen.
Wie das implizite Euler-Verfahren funktioniert weiß ich. Da ich aber 4 Gleichungen habe, die zusammenhängen, kann ich nicht einfach nach dv(i) auslösen. Kann mir jemand sagen wie ich dieses Problem lösen kann?

Liebe Grüße,
Mieze

Hier mein Code:

Code:

close all
clear all
clc

% Anfangswerte / Parameter
dV_0 = -65;
betta = 0;
t = 80;
N = 8000;

% zusätzlich für HH
C = 1;
g_Na = 120;
g_K = 36;
g_L = 0.3;
E_Na = 50;
E_K = -77;
E_L = -54.4;

% Funktion aufrufen
[dV,dn,dm,dh] = HodgkinHuxley77(dV_0, betta, N, t, C, g_Na, g_K, g_L, ...
E_Na, E_K, E_L, I_app);

function [ dV,dn,dm,dh ] = HodgkinHuxley77( dV_0, betta, N, t, C, g_Na, g_K, g_L, E_Na, E_K, E_L, I_app)

delta_t = t/N; %Schrittweite

%Hodgkin-Huxley coefficients
alpha_n = @(dV) 0.01*(dV + 55)/(1-exp(-(dV + 55)/10));
beta_n = @(dV) 0.125*exp(-(dV + 65)/80);
alpha_m = @(dV) 0.1*(dV + 40)/(1 - exp(-(dV +40)/10));
beta_m = @(dV) 4*exp(-(dV + 65)/18);
alpha_h = @(dV) 0.07*exp(-(dV + 65)/20);
beta_h = @(dV) 1/(1 + exp(-(dV + 35)/10));

n_inf = @(dV) alpha_n(dV)/(alpha_n(dV) + beta_n(dV));
tau_n = @(dV) 1/(alpha_n(dV) + beta_n(dV));
m_inf = @(dV) alpha_m(dV)/(alpha_m(dV) + beta_m(dV));
tau_m = @(dV) 1/(alpha_m(dV) + beta_m(dV));
h_inf = @(dV) alpha_h(dV)/(alpha_h(dV) + beta_h(dV));
tau_h = @(dV) 1/(alpha_h(dV) + beta_h(dV));

dV = zeros(1,N+1); % V als Matrix
dn = zeros(1,N+1); % n als Matrix
dm = zeros(1,N+1); % m als Matrix
dh = zeros(1,N+1); % h als Matrix

%Anfangswerte festlegen
dV(:,1) = dV_0;
dn(:,1) = alpha_n(dV_0)/(alpha_n(dV_0) + beta_n(dV_0));%dn_0;
dm(:,1) = alpha_m(dV_0)/(alpha_m(dV_0) + beta_m(dV_0));%dm_0;
dh(:,1) = alpha_h(dV_0)/(alpha_h(dV_0) + beta_h(dV_0));%dh_0;

Z = randn(100,N+1); % entspricht N_k, Zufallsvariable

for i = 2:N+1
for j = 1
dV(j,i) = dV(j,i-1) + ((1/C)*((I_app(i) - g_K*(dn(j, i))^4*(dV(j,i) - E_K) - g_Na*(dm(j, i))^3*dh(j, i)*(dV(j,i) - E_Na) - g_L*(dV(j,i) - E_L) + betta*sqrt(delta_t)*Z(j,i))*delta_t ));
dn(j,i) = dn(j,i-1) + (((n_inf(dV(j,i)) - dn(j, i))/tau_n(dV(j, i)))*delta_t);
dm(j,i) = dm(j,i-1) + (((m_inf(dV(j,i)) - dm(j, i))/tau_m(dV(j, i)))*delta_t);
dh(j,i) = dh(j,i-1) + (((h_inf(dV(j,i)) - dh(j, i))/tau_h(dV(j, i)))*delta_t);
end
end

Funktion ohne Link?

PS: Meine Frage bezieht sich auf die 4 Gleichungen in der for-Schleife, also dV, dn, dm und dh.

Harald

Forum-Meister


	Beiträge: 24.499

	Anmeldedatum: 26.03.09

	Wohnort: Nähe München

	Version: ab 2017b

Verfasst am: 19.06.2015, 21:26 Titel:

Hallo,

implizite Verfahren beinhalten, dass in jedem Iterationsschritt ein (generell) nichtlineares Gleichungssystem gelöst in $x_{k+1}$ werden muss:
$x_k+h * f(t_{k+1},x_{k+1}) - x_{k+1} = 0;$

Das kann z.B. mit fsolve gemacht werden. Die Dimensionen von $x_{k+1}$ sind dafür recht gleichgültig.

Der Wer
Grüße,
Harald

mieze579

Themenstarter

Forum-Newbie


	Beiträge: 6

	Anmeldedatum: 22.01.15

	Wohnort: ---

	Version: ---

Verfasst am: 20.06.2015, 13:33 Titel:

Danke für die schnelle Antwort. Ich kann die Gleichungen

Code:

for i = 2:N+1
for j = 1
dV(j,i) = dV(j,i-1) + ((1/C)*((I_app(i) - g_K*(dn(j, i))^4*(dV(j,i) - E_K) - g_Na*(dm(j, i))^3*dh(j, i)*(dV(j,i) - E_Na) - g_L*(dV(j,i) - E_L) + betta*sqrt(delta_t)*Z(j,i))*delta_t ));
dn(j,i) = dn(j,i-1) + (((n_inf(dV(j,i)) - dn(j, i))/tau_n(dV(j, i)))*delta_t);
dm(j,i) = dm(j,i-1) + (((m_inf(dV(j,i)) - dm(j, i))/tau_m(dV(j, i)))*delta_t);
dh(j,i) = dh(j,i-1) + (((h_inf(dV(j,i)) - dh(j, i))/tau_h(dV(j, i)))*delta_t);
end
end

Funktion ohne Link?

ja so umschreiben, dass auf der linken Seite eine Null steht, also

Code:

WertV=0;
Wertn=0;
Wertm=0;
Werth=0;

for i = 2:N+1
for j = 1
WertV= -dV(j,i)+ dV(j,i-1) + ((1/C)*((I_app(i) - g_K*(dn(j, i))^4*(dV(j,i) - E_K) - g_Na*(dm(j, i))^3*dh(j, i)*(dV(j,i) - E_Na) - g_L*(dV(j,i) - E_L) + betta*sqrt(delta_t)*Z(j,i))*delta_t ));
Wertn = -dn(j,i)+dn(j,i-1) + (((n_inf(dV(j,i)) - dn(j, i))/tau_n(dV(j, i)))*delta_t);
Wertm = -dm(j,i) + dm(j,i-1) + (((m_inf(dV(j,i)) - dm(j, i))/tau_m(dV(j, i)))*delta_t);
Werth = -dh(j,i) + dh(j,i-1) + (((h_inf(dV(j,i)) - dh(j, i))/tau_h(dV(j, i)))*delta_t);
end
end

Funktion ohne Link?

Und wie genau wende ich jetzt fsolve an?Für eine "einfache" Gleichung würde ich die Anwendung verstehen, allerdings weiß ich nun nicht, wie ich die Anfangswerte für dV, dn, dm und dh einfließen lasse.
Brauche ich denn pro Gleichung einmal fsolve oder muss ich eine größere Matrix mit allen 4 Gleichungen bauen?

Harald

Forum-Meister


	Beiträge: 24.499

	Anmeldedatum: 26.03.09

	Wohnort: Nähe München

	Version: ab 2017b

Verfasst am: 20.06.2015, 13:46 Titel:

Hallo,

Zitat:

Brauche ich denn pro Gleichung einmal fsolve oder muss ich eine größere Matrix mit allen 4 Gleichungen bauen?

Letzteres, denn die Gleichungen hängen ja zusammen.

Wenn du die rechte Seite der DGL erstmal in eine Funktion f auslagerst, dürfte es klarer werden. Dann hast du

Code:

x(:, k+1) = fsolve(@(x) x(:,k) + f(t(k+1), x) - x, x(:,k);

Funktion ohne Link?

Grüße,
Harald

mieze579

Themenstarter

Forum-Newbie


	Beiträge: 6

	Anmeldedatum: 22.01.15

	Wohnort: ---

	Version: ---

Verfasst am: 23.06.2015, 14:23 Titel:

Ok, aber so ganz habe ich das noch nicht verstanden.

Deine Berechnung:

Code:

x(:, k+1) = fsolve(@(x) x(:,k) + f(t(k+1), x) - x, x(:,k);

Funktion ohne Link?

.
In deiner Berechnung fehlt eine Klammer. Und wieso nehme ich nur f(t(k+1), x) und nicht f(t(k+1), x+1)? So wäre es doch im impliziten Euler-Verfahren. Und wieso wird danach noch ein x abgezogen?
Das x(:,k) am Ende ist dann der Anfangswert?Müsste ich da nicht einen Vektor einfügen, der die Anfangswerte für V, n, m und h enthält?

Ich habe nun die Berechnung von V, n, m und h in eine Matrix zusammengefügt und in eine Funktion ausgelagert:

Code:

function Matrix = myfun(N, dV, dn, dm, dh, C, I_app, g_K, g_L, g_Na, E_K, E_L,E_Na, betta, delta_t, Z, n_inf, m_inf, h_inf, tau_n, tau_m, tau_h)
Matrix = zeros(N+1,4);
Matrix(1, :) = [dV(:,1), dn(:,1), dm(:,1), dh(:,1)];
for i = 2:N+1
%beginne in erster Zeite! M=(dV, dn, dm, dh) (jewileige Vektoren in Spalten von M)
Matrix(i-1, :) = [(-dV(i) + dV(i-1) + ((1/C)*((I_app(i) - g_K*(dn(i))^4*(dV(i) - E_K) - g_Na*(dm(i))^3*dh(i)*(dV(i) - E_Na) - g_L*(dV(i) - E_L) + betta*sqrt(delta_t)*Z(i))*delta_t )));
-dn(i) + dn(i-1) + (((n_inf(dV(i)) - dn(i))/tau_n(dV(i)))*delta_t);
-dm(i) + dm(i-1) + (((m_inf(dV(i)) - dm(i))/tau_m(dV(i)))*delta_t);
-dh(i) + dh(i-1) + (((h_inf(dV(i)) - dh(i))/tau_h(dV(i)))*delta_t)];
end

Funktion ohne Link?

Nun sind in Matrix die Vektoren V, n, m und h jeweils in den Spalten angeordnet.
Jetzt rufe ich diese Funktion wieder auf:

Code:

Matrix2 = zeros(N+1, 4);
for k = 1:N
Matrix2(:, k+1) = fsolve(@myfun, Matrix(:,k));
end

Funktion ohne Link?

Wenn ich jetzt alles durchlaufen lasse kommen einige Fehler:
Warning: Input arguments must be scalar.
> In myfun at 7
In fsolve at 254
In HodgkinHuxley77 at 56
In Aufruf77 at 55
??? Input argument "dV" is undefined.

Error in ==> myfun at 8
Matrix(1, Smile

= [dV(:,1), dn(:,1), dm(:,1), dh(:,1)];

Error in ==> fsolve at 254
fuser = feval(funfcn{3},x,varargin{:});

Error in ==> HodgkinHuxley77 at 56
Matrix2(:, k+1) = fsolve(@myfun, Matrix(:,k));

Error in ==> Aufruf77 at 55
[dV,dn,dm,dh] = HodgkinHuxley77(dV_0, betta, N, t, C,
g_Na, g_K, g_L, ...

Caused by:
Failure in initial user-supplied objective function
evaluation. FSOLVE cannot continue.

Kannst du mir da helfen?

Harald

Forum-Meister


	Beiträge: 24.499

	Anmeldedatum: 26.03.09

	Wohnort: Nähe München

	Version: ab 2017b

Verfasst am: 23.06.2015, 18:52 Titel:

Hallo,

Zitat:

In deiner Berechnung fehlt eine Klammer.

Stimmt; ganz am Ende.

Zitat:

Und wieso nehme ich nur f(t(k+1), x) und nicht f(t(k+1), x+1)? So wäre es doch im impliziten Euler-Verfahren.

Mhh, das war unglücklich. So sollte es verständlicher und richtiger sein:

Code:

x(:, k+1) = fsolve(@(xNext) x(:,k) + f(t(k+1), xNext) - xNext, x(:,k));

Funktion ohne Link?

xNext ist hier eine Art Dummy-Variable. Die Lösung wird dann in die (k+1)-te Spalte von x geschrieben.

Zitat:

Und wieso wird danach noch ein x abgezogen?

Implizite Euler-Formel, und alles auf eine Seite gebracht. Siehe auch Beitrag 19.06.2015, 21:26.

Zur Nutzung von fsolve an sich: in der Doku zu fsolve ist folgendes verlinkt, bitte mal insbesondere den Teil zu "Anonymous Functions" durcharbeiten:
http://de.mathworks.com/help/optim/.....ing-extra-parameters.html

Ansonsten blicke ich bei der Struktur der Funktion noch nicht ganz durch.
Das f aus dem obigen Code sollte die rechte Seite der DGL sein:
y' = f(t, y)

Dabei kann y durchaus ein Vektor sein, es sollte aber keine Matrix sein.

Grüße,
Harald


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 - 2025 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.