ich versuche gerade, eine eindimensionale Fokker.Planck Gleichung durch Finite Differenzen zu lösen. die Fokker-Planck Gleichung beschreibt die Ausbreitung der Dichte eines stochastischen Prozesses. ich muss eigentlich nur diese Gleichung durch Diskretisierung lösen. Zuerst in x-Richtung, danach in t-Richtung.
Als Randbedingung soll ich eine Dirichlet- Bedingung oder eine Robin Bedingung versuchen.
Als Anfangsbedingung zum Zeitpunkt t=0 soll eine Standardnormalverteilung zugrunde liegen.
Zusätzlich soll f natürlich eine Dichte sein und somit
. Dazu weiß ich nicht, wie ich das MATLAB beibringe..
Hier ist mein Code
Ich hab bis jetzt folgendes versucht.
Code:
%Finite Differenzen zur Fokker- Planck Gleichung
%clear; clc;
%Initialisierung
nx=60; %Anzahl der Gitterpunkte in x-Richtung.
nt=12; %Anzahl der Gitterpunkte in t-Richtung.
T=1; %Endzeitpunkt.
X= 10; %Länge der x-Achse.
h=1/(nx-1); %Schrittweite x-Richtung.
dt=T*h^2; %Schrittweite t-Richtung.
sigma=0.12;
%Anfangsbedingung, f Vektor mit Standard-Normalverteilungen
x0=0;
st=1;
%x = [0:h:100] for j=1:nx
f = ones(nx,1) %*(1/sqrt(2*pi)*st)*exp(-0.5*((j*h-x0)/st)^2);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Differenzenapproximation. sowohl mu als auch f werden diskretisiert. % Zuerst nach dem Ort, dann nach der Zeit.
t = 0;
%Versuch eines expliziten Verfahrens
for k=1:nt
falt=f;
t=t+dt;
f(1)=0;
mu(nx)=-0.015*(nx*h-2);
f(nx)=0;
for i=2:nx-1 %für jede Zeitscheibe, löse FKP nach f(i) unter f(i-1) zur Zeit t auf
mu(i)=-0.015*(i*h-2);
f(i)= falt(i) + dt*((1/2)*h*(-mu(i+1)*falt(i+1)-mu(i-1)*falt(i-1))+ 0.5*sigma^2*(falt(i+1)-2*falt(i)+falt(i-1))*(1/h^2)) end
An der Steller %%%%%%%%%%%%%%%%%%%%%% weiß ich nicht, wie ich die Anfangsbedingung, dass f eine Funktion, also die Standardnormalverteilung von x ist, integrieren kann.ich muss sie ja auch, um die Diskretisierungsschleife durchzuführen, als (nx,1) Vektor definieren. Ich möchte, dass f0=(1/sqrt(2*pi)*1)*exp(-0.5*((x-0)/1)^2) ist und dann bei der Diskretisierung immer an den jeweiligen Stellen x=i*h ausgewertet wird..Leider kein Plan wie..
Ich weiß hier fehlt noch einiges, aber ich wäre echt sau dankbar, wenn mir jemand helfen könnte... Bin leider nen MATLAB Anfänger, wie ihr seht...
Zuletzt bearbeitet von Fettknobelchen85 am 27.11.2011, 14:55, insgesamt einmal bearbeitet
Hey! Ich habe noch eine zweite Möglichkeit problert, die etwas vielversprechender aussieht. Vielleicht wollt ihr euch das auch noch ansehen, die Probleme sind die gleichen...
Code:
%Neuer Versuch Finite Differenzen
%clear; clc;
X = 50; %Länge in x Richtung
T = 80; %Endzeitpunkt
tx = 50;
dt = T/tx;
nx = 60;
sigma= 0.12;
h = X/nx;
b = dt/(h^2);
% % Stützstellen x=i*h % for i=1:n-1 % x(i)=(i+1)*h; % end
Wie gesagt, ich wäre unendlich dankbar, wenn sich das mal jemand ansehen könnte und mir tipps geben könnte, wie ich die Anfangs und Randbedingungen und, die integral über fkt =1 Bedingung einfüge. Hab das hier in der for schleife versucht über u0. hier hab ich es auch mit 2 input argumente der Fukntion, (die ich hier u genannt hab, damit ich nicht durcheinander komme mit den Programmen) versucht.
Soll insgesamt ein explizites Finite Differenzen Verfahren geben. Hier noch einmal die Ausgangsgleichung, eine eindimensionale parabolische partielle Differentialgleichung 2. Ordnung.
Das sigma sei konstant angenommen und mu eine lineare Funktion, hier -0.015*(x-2).
Falls du die Pdgl nicht per hand lösen willst würde ich die Funktion pdepe empfehlen.
Da die dgl kein Quellterm hat sollte doch das integral über f konstant sein (falls die RB entsprechend gewählt sind). Diese Bedingung ist also automatisch erfüllt.
Danke für die Antwort. Mit PDEPE habe ich es auch schon versucht. Das Problem ist, dass ich aus der PDE ein System von ODEs erstellen muss. Dazu diskretisiere ich das von Hand.
Also ich muss es leider so machen... Aus den Daten, die ich dabei erzeuge will ich nachher rückwärts die Parameter der Funktion mu schätzen...
Im Code fehlt übrigens noch mu=zeros(nx,1);
Das das Integral konstant bleibt dachte ich mir auch schon, dazu muss ich aber die Anfangbedingung integrieren. Also das
u0= (1/sqrt(2*pi)*st)*exp(-0.5*((x-x0)/st)^2)
das ist in der schleife glaub ich noch falsch.. Wie bekomme ich das denn hin, dass ich dort ein x verwenden kann? ich muss MATLAB doch sagen, dass dieses x an den Diskretisierungsstellen ausgewertet werden soll...
Wie ich die Randbedingungen nun aufstelle, weiß ich leider auch noch nicht... Es würde schon eine einfache Dirichlet Bedingung reichen..
hier mal ein paar verbesserungen. funktioniert aber noch nicht ganz.
Um die Anfangsbedingungen zu implementieren brauchst du keine Schleife, da man in Matlab sehr elegant mit Matrizen rechnen kann.
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
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.
RSS