A simple example of coding and solving the Fokker-Planck equation with the Finite Element Method in JAVA

Suzuki Takashi (October 20, 2010 12:23 PM) | コメント(0)

At first I will be touching on the correspondence of the Fokker-Planck equation to the stochastic differential equation.
$$ dX_t\ =\ D_1dt+D_2dW_t$$

The above stochastic differential equation corresponds to the below Fokker-Planck equation.
$$ \frac{\partial}{\partial t}f(x,t)=-\frac{\partial}{\partial x}D_1f(x,t)+\frac{1}{2}\frac{\partial^2}{\partial x^2}D_2^2f(x,t)$$

We will be solving the above Fokker-Planck equation analytically using the below transformation into the below easy problem.
$$ f(x,t)=\exp\left\{\frac{D_1\left(x-\frac{D_1t}{2}\right)}{D_2^2}\right\} g(x,t)$$
$$ \frac{\partial}{\partial t}g(x,t)=\frac{1}{2}\frac{\partial^2}{\partial x^2}D_2^2g(x,t)$$
$$ g(x,t)=\frac{1}{\sqrt{2\pi D_2^2 t}}\exp\left(-\frac{x^2}{2D_2^2t}\right)$$

We will acquire the below solution to the above Fokker-Planck equation as a result of the above transformation.
$$ f(x,t)=\frac{1}{\sqrt{2\pi D_2^2 t}}\exp\left(-\frac{(x-D_1t)^2}{2D_2^2t}\right)$$

We will find the theoretical value of the above Fokker-Planck equation using the above solution.

Well, I will be solving the above Fokker-Planck equation with the finite element method. Let the tent function and f(x,t) be defined as the below equation.
$$\phi_i(x)= \begin{cases} \displaystyle \frac {x-x_{i-1}}{h} & \left(x_{i-1}<x\leqq x_i\right) \\ \displaystyle \frac{x_{i+1}-x}{h} & \left(x_i<x\leqq x_{i+1}\right) \end{cases}$$
$$f(x,t)= \begin{cases} f_i\phi_i(x)+f_{i-1}\phi_{i-1}(x) & \left(x_{i-1}<x\leqq x_i\right) \\ f_i\phi_i(x)+f_{i+1}\phi_{i+1}(x) & \left(x_i<x\leqq x_{i+1}\right) \end{cases}$$

I will also describe the Fokker-Planck equation to use the finite element method.
$$ \int_{x_0}^{x_{n+1}}\left(\frac{\partial f}{\partial t}+D_1\frac{\partial f}{\partial x}-\frac{1}{2}D_2^2\frac{\partial^2f}{\partial x^2}\right)\phi_idx=0$$
$$ \int_{x_0}^{x_{n+1}}\dfrac{\partial f}{\partial t}\phi_idx+D_1\int_{x_0}^{x_{n+1}}\dfrac{\partial f}{\partial x}\phi_idx+\dfrac{1}{2}D_2^2\int_{x_0}^{x_{n+1}}\dfrac{\partial f}{\partial x}\dfrac{\partial \phi_i}{\partial x}dx=0 $$
$$\begin{aligned}&\left\{\dfrac{h}{6\Delta t}+\left(-\dfrac{D_1}{2}-\dfrac{D_2^2}{2h}\right)\right\}f_{i-1}^{k+1}+\left(\dfrac{4h}{6\Delta t}+\dfrac{D_2^2}{h}\right)f_i^{k+1} \\ &+\left\{\dfrac{h}{6\Delta t}+\left(\dfrac{D_1}{2}-\dfrac{D_2^2}{2h}\right)\right\}f_{i+1}^{k+1} = \dfrac{h}{6\Delta t}f_{i-1}^k+\dfrac{4h}{6\Delta t}f_i^k+\dfrac{h}{6\Delta t}f_{i+1}^k\end{aligned}$$

I will represent the above Fokker-Planck equation as the below matrix form.
$$ Uf^{k+1}\ =\ Vf^k$$
$$ f^{k+1}=\begin{pmatrix} f_1^{k+1} & f_2^{k+1} & f_3^{k+1} & \cdots & f_n^{k+1} \end{pmatrix}^T$$
$$ f^{k}=\begin{pmatrix} f_1^{k} & f_2^{k} & f_3^{k} & \cdots & f_n^{k} \end{pmatrix}^T$$
$$U=\begin{pmatrix} \dfrac{4h}{6\Delta t}+\dfrac{D_2^2}{h} & \dfrac{h}{6\Delta t}+\left(\dfrac{D_1}{2}-\dfrac{D_2^2}{2h}\right) & 0 & \cdots & 0 \\[1.5ex] \dfrac{h}{6\Delta t}+\left(-\dfrac{D_1}{2}-\dfrac{D_2^2}{2h}\right) & \dfrac{4h}{6\Delta t}+\dfrac{D_2^2}{h} & \dfrac{h}{6\Delta t}+\left(\dfrac{D_1}{2}-\dfrac{D_2^2}{2h}\right) & \cdots & 0 \\[1.5ex] 0 & \dfrac{h}{6\Delta t}+\left(-\dfrac{D_1}{2}-\dfrac{D_2^2}{2h}\right) & \dfrac{4h}{6\Delta t}+\dfrac{D_2^2}{h} & \cdots & 0 \\[1.5ex] \vdots & \vdots & \vdots & \ddots & \vdots \\[1.5ex] 0 & 0 & 0 & \cdots & \dfrac{4h}{6\Delta t}+\dfrac{D_2^2}{h} \end{pmatrix}$$
$$V=\begin{pmatrix} \dfrac{4h}{6\Delta t} & \dfrac{h}{6\Delta t} & 0 & \cdots & 0 \\[1.5ex] \dfrac{h}{6\Delta t} & \dfrac{4h}{6\Delta t} & \dfrac{h}{6\Delta t} & \cdots & 0 \\[1.5ex] 0 & \dfrac{h}{6\Delta t} & \dfrac{4h}{6\Delta t} & \cdots & 0 \\[1.5ex] \vdots & \vdots & \vdots & \ddots & \vdots \\[1.5ex] 0 & 0 & 0 & \cdots & \dfrac{4h}{6\Delta t} \end{pmatrix}$$

U and V are the tridiagonal matrices. The example of the finite element method is given below.

      source code: FokkerPlanck.java

import java.awt.*;
import java.awt.event.*;
import java.text.DecimalFormat;
import java.applet.Applet;
import javax.swing.JLabel;
import javax.swing.JTextField;

public class FokkerPlanck extends Applet implements ActionListener{

	// variable set
	private static final long serialVersionUID = 9129736202608299398L;
	public static double d1 = 2;
	public static double d2 = 3;
	public static double t = 1;
	public static double h = 0.055;
	public static double dt = 0.25;
	public static double[] s1 = new double[1000];
	public static double[] s2 = new double[1000];
	public static double[] s3 = new double[1000];
	public static double[] s4 = new double[1000];
	public static int[][] xxxx = new int [56][1002];
	public static int[][] yyyy = new int [56][1002];
	public static int[] xxxx1 = new int [1002];
	public static int[] yyyy1 = new int [1002];
	public static int[] xxxx2 = new int [1000];
	public static int[] yyyy2 = new int [1000];
	public static int[] xxxx3 = new int [1000];
	public static int[] yyyy3 = new int [1000];
	public static double[][] MatrixL = new double[1000][1000];
	public static double[][] MatrixU = new double[1000][1000];
	public static double[][] Matrix01 = new double[1000][1000];
	public static double[][] Matrix02 = new double[1000][1000];
	public static double[][] Matrix03 = new double[1000][1000];
	public static double[][] vector01 = new double[57][1000];
	public static double[][] vector02 = new double[56][1000];
	public static double[][] vector03 = new double[56][1000];
	JTextField yx = new JTextField("2.0");
	JTextField yn = new JTextField("3.0");
	JTextField yr = new JTextField("0.055");
	JLabel label1 = new JLabel("The Example of The Finite Element Method",JLabel.CENTER);
	JLabel label2 = new JLabel("<html><body><font size=3>The Drift Coefficient (0<D<sub>1</sub>≤10<sup>8</sup>)</font></body></html>",JLabel.CENTER);
	JLabel label3 = new JLabel("<html><body><font size=3>The Diffusion Coefficient (0<D<sub>2</sub>≤10<sup>8</sup>)</font></body></html>",JLabel.CENTER);
	JLabel label4 = new JLabel("<html><body><font size=3>The Interval between x<sub>i</sub> and x<sub>i+1</sub></font></body></html>", JLabel.CENTER);

	public void init(){
		label1.setPreferredSize(new Dimension(416,24));
		label1.setFont(new Font("Serif",Font.BOLD,14));
		add(label1);
		label2.setPreferredSize(new Dimension(260,23));
		label2.setFont(new Font("Serif",Font.BOLD,11));
		add(label2);
		yx.setPreferredSize(new Dimension(80,23));
		add(yx);
		label3.setPreferredSize(new Dimension(260,23));
		label3.setFont(new Font("Serif",Font.BOLD,11));
		add(label3);
		yn.setPreferredSize(new Dimension(80,23));
		add(yn);
		label4.setPreferredSize(new Dimension(260,23));
		label4.setFont(new Font("Serif",Font.BOLD,11));
		add(label4);
		yr.setPreferredSize(new Dimension(80,23));
		add(yr);
		yx.addActionListener(this);
		yn.addActionListener(this);
		yr.addActionListener(this);
	}

	public void actionPerformed(ActionEvent e){
		if(e.getSource()==yx){
			d1=Double.valueOf(yx.getText()).doubleValue();
			if(d1 <= 0 || d1 > Math.pow(10.0,8)){
				d1 = 2.0;
			}
		}
		if(e.getSource()==yn){
			d2 = Double.valueOf(yn.getText()).
			doubleValue();
			if(d2 <= 0 || d2 > Math.pow(10.0,8)){
				d2 = 3.0;
			}
		}
		if(e.getSource()==yr){h = Double.valueOf
		(yr.getText()).doubleValue();
			if(h <= 0 || h > Math.pow(10.0,8)){h = 0.055;}
		}
		yx.setText(""+d1);
		yn.setText(""+d2);
		yr.setText(""+h);
		repaint();
	}

	public void paint(Graphics g){

		for(int i=0;i<1000;i++){
			s1[i]=(i*0.05-10)*d2/3;
		}
		for(int i=0;i<1000;i++){
			s2[i]=function(d1,d2,s1[i],1);
		}
		for(int i=0;i<1000;i++){
			s3[i]=function(d1,d2,s1[i],6);
		}
		for(int i=0;i<1000;i++){
			s4[i]=function(d1,d2,s1[i],14);
		}

		Matrix01[0][0]=4*h/(6*dt)+Math.pow(d2,2)/h;
		Matrix01[0][1]=h/(6*dt)+d1/2-Math.pow(d2,2)/(2*h);
		for(int i=1;i<999;i++){
			Matrix01[i][i-1]=h/(6*dt)-d1/2-Math.pow(d2,2)/(2*h);
			Matrix01[i][i]=4*h/(6*dt)+Math.pow(d2,2)/h;
			Matrix01[i][i+1]=h/(6*dt)+d1/2-Math.pow(d2,2)/(2*h);
		}
		Matrix01[999][998]=h/(6*dt)-d1/2-Math.pow(d2,2)/(2*h);
		Matrix01[999][999]=4*h/(6*dt)+Math.pow(d2,2)/h;

		Matrix02[0][0]=4*h/(6*dt);
		Matrix02[0][1]=h/(6*dt);
		for(int i=1;i<999;i++){
			Matrix02[i][i-1]=h/(6*dt);
			Matrix02[i][i]=4*h/(6*dt);
			Matrix02[i][i+1]=h/(6*dt);
		}
		Matrix02[999][998]=h/(6*dt);
		Matrix02[999][999]=4*h/(6*dt);

		vector01[0][0]=Matrix02[0][0]*s2[0]+Matrix02[0][1]*s2[1];
		for(int i=1;i<999;i++){
			vector01[0][i]=Matrix02[i][i-1]*s2[i-1]+Matrix02[i][i]*s2[i]+Matrix02[i][i+1]*s2[i+1];
		}
		vector01[0][999]=Matrix02[999][998]*s2[998]+Matrix02[999][999]*s2[999];

                for(int i=0;i<1000;i++){
                    for(int j=0;j<1000;j++){
                        MatrixU[i][j]=Matrix01[i][j];
                    }
                }
                for(int i=1;i<1000;i++){
                    MatrixU[i][i]=Matrix01[i][i]-MatrixU[i-1][i]*Matrix01[i][i-1]/MatrixU[i-1][i-1];
                }
                MatrixL[0][0]=1;
                for(int i=1;i<1000;i++){
                    MatrixL[i][i-1]=Matrix01[i][i-1]/MatrixU[i-1][i-1];
                    MatrixL[i][i]=1;
                }

		for(int k=0;k<56;k++){
			vector02[k][0]=vector01[k][0];
			for(int i=1;i<1000;i++){
				vector02[k][i]=vector01[k][i]-MatrixL[i][i-1]*vector02[k][i-1];
			}
			vector03[k][999]=vector02[k][999]/MatrixU[999][999];
			for(int i=998;i>=0;i--){
				vector03[k][i]=(vector02[k][i]-MatrixU[i][i+1]*vector03[k][i+1])/MatrixU[i][i];
			}
			vector01[k+1][0]=Matrix02[0][0]*vector03[k][0]+Matrix02[0][1]*vector03[k][1];
			for(int i=1;i<999;i++){
				vector01[k+1][i]=Matrix02[i][i-1]*vector03[k][i-1]+Matrix02[i][i]*vector03[k][i]+Matrix02[i][i+1]*vector03[k][i+1];
			}
			vector01[k+1][999]=Matrix02[999][998]*vector03[k][998]+Matrix02[999][999]*vector03[k][999];
		}

		double	minss1 = s1[0];
		double	maxss1 = s1[0];

		for(int i=0;i<1000;i++){
			if(s1[i]>maxss1){
				maxss1 = s1[i];
			}
		}
		for(int i=0;i<1000;i++){
			if(s1[i]<minss1){
				minss1 = s1[i];
			}
		}

		for(int j=0;j<56;j++){
			for(int i=0;i<1000;i++){
				xxxx[j][i] = (int)(i*351/999)+33;
				yyyy[j][i] = 295-(int)((vector03[j][i])*176/0.2);
			}
			xxxx[j][1000]=33;
			xxxx[j][1001]=387;
			yyyy[j][1000]=299;
			yyyy[j][1001]=299;
		}

		for(int i=0;i<1000;i++){
			xxxx1[i] = (int)(i*351/999)+33;
			yyyy1[i] = 295-(int)((s2[i])*176/0.2);
		}
		xxxx1[1000]=33;
		xxxx1[1001]=387;
		yyyy1[1000]=299;
		yyyy1[1001]=299;

                for(int i=0;i<1000;i++){
			xxxx2[i] = (int)(i*351/999)+33;
			yyyy2[i] = 295-(int)((s3[i])*176/0.2);
		}

		for(int i=0;i<1000;i++){
			xxxx3[i] = (int)(i*351/999)+33;
			yyyy3[i] = 295-(int)((s4[i])*176/0.2);
		}

		Graphics2D g2 = (Graphics2D)g;
		GradientPaint gp1 = new GradientPaint(0, 0, new Color(154,181,228), 0,470,new Color(225,232,245), true);
		g2.setPaint(gp1);
		g2.fillRect(0,0,416,328);
		super.paint(g);
		GradientPaint gp2 = new GradientPaint(0, 33, new Color(225,232,245), 0,351,new Color(154,181,228), true);
		g2.setPaint(gp2);
		g2.fillRect(33,123,354,176);
		GradientPaint gp4 = new GradientPaint(33,227, new Color(192,88,77), 33,418,new Color(160,82,45), true);
		g2.setPaint(gp4);
		for(int k=0;k<2;k++){
			for(int i=0;i<12;i++){
				g2.setPaint(gp4);
				g2.fillPolygon(xxxx1,yyyy1,1001);
			}
			g2.setPaint(gp2);
			g2.fillRect(33,123,354,176);
			for(int i=0;i<56;i+=4){
				for(int j=0;j<12;j++){
				g2.setPaint(gp4);
				g2.fillPolygon(xxxx[i],yyyy[i],1001);
				}
				g2.setPaint(gp2);
				g2.fillRect(33,123,354,176);
                        }
			g2.setPaint(gp4);
			g2.fillPolygon(xxxx[55],yyyy[55],1001);
			g2.setPaint(gp2);
			g2.fillRect(33,123,354,176);
		}
		for(int i=0;i<1000;i++){
				g2.setColor (new Color(0,0,0));
				g2.fillOval(xxxx1[i], yyyy1[i], 3, 3);
		}
		for(int i=0;i<1000;i++){
				g2.fillOval(xxxx2[i], yyyy2[i], 3, 3);
		}
		for(int i=0;i<1000;i++){
				g2.fillOval(xxxx3[i], yyyy3[i], 3, 3);
		}
		for(int j=0;j<1000;j++){
				g2.setPaint(gp4);
				g2.fillOval(xxxx[23][j], yyyy[23][j], 3, 3);
				g2.fillOval(xxxx[55][j], yyyy[55][j], 3, 3);
		}
		GradientPaint gp3 = new GradientPaint(33, 227, new Color(79,129,189), 33,418,new Color(65,105, 225), true);
		g2.setPaint(gp3);
		DecimalFormat df0 = new DecimalFormat("0.00");
		DecimalFormat df1 = new DecimalFormat("0.000");
		for(int i=0;i<5;i++){
			g2.setFont(new Font ("Serif",Font.BOLD,12));
			g2.setColor (new Color(0,0,0));
			String subx = df1.format((maxss1-minss1)/4*i+minss1);
			g2.drawString(subx,14+88*i,313);
		}
		for(int i=0;i<3;i++){
			g2.setFont(new Font ("Serif",Font.BOLD,12));
			g2.setColor (new Color(0,0,0));
			String suby = df0.format(0.1*i);
			g2.drawString(suby,3,299-85*i);
		}
		g2.setFont(new Font ("Serif",Font.BOLD,12));
		g2.setColor (new Color(0,0,0));
		g2.drawString("n = 1000",250,152);
		g2.setColor (new Color(79,129,189));
		BasicStroke BoldLine = new BasicStroke(2.8f);
		g2.setStroke(BoldLine);
		g2.drawLine(250,182,320,162);
		g2.setColor (new Color(0,0,0));
		g2.drawString("The Analytic Solution",250,182);
		g2.setPaint(gp4);
		g2.drawLine(250,212,320,192);
		g2.setColor (new Color(0,0,0));
		g2.drawString("The FEM Solution",250,212);
		g2.drawLine(250,222,320,222);
		g2.setColor (new Color(0,0,0));
		g2.drawString("The Initial Value",250,242);
		g2.drawString("f",15,115);
		g2.drawString("x",340,313);
	}

	public static double function(double d1, double d2, double x, double t){
		return 1/Math.pow(Math.PI*2*Math.pow(d2,2)*t,0.5)*Math.exp(-Math.pow((x-d1*t),2)/(2*Math.pow(d2,2)*t));
	}
}

In conclusion the difference between the analytic solution and the solution with the finite element method seems to depend on the form of the tent function and other parameters and still more improvement seems to be necessary. There is a great deal of room for improvement.

Finally I am happy to assist you in coding and solving the Fokker-Planck equation with the finite element method in JAVA.

カテゴリ:

コメントする

My Photo
プロフィール！
2016･11･15 改訂
spacer01
rssspacer01foaf
spacer01
atom.xml
spacer01

この記事について

このページは、Suzuki TakashiがOctober 20, 2010 12:23 PMに書いた記事です。

ひとつ前の記事は「調査再開、パート8。」です。

次の記事は「まだ調査が必要かな。」です。

最近のコンテンツはインデックスページで見られます。過去に書かれたものはアーカイブのページで見られます。

June 2024

Sun Mon Tue Wed Thu Fri Sat
1
2 3 4 5 6 7 8
9 10 11 12 13 14 15
16 17 18 19 20 21 22
23 24 25 26 27 28 29
30

カテゴリ

月別 アーカイブ

OpenID対応しています OpenIDについて