泊松-费米方程数值求解

常见于离子液体等带电软物质体系的泊松-费米方程

\begin{equation} \frac{d^2\phi}{dx^2}=\frac{\sinh \phi}{1+2\gamma \sinh^2 \phi/2}=\rho(\phi) \label{Poisson-Fermi2} \end{equation}

边界条件：

\begin{equation} \begin{split} \phi(0)=&V_0 \\ \phi(\infty)=&0 \end{split} \label{BC2} \end{equation}

下面我们用bvp4c 解方程，重复出文献 Double Layer in Ionic Liquids: Overscreening versus Crowding中 FIG. 2(a) 中的虚线。

依次取 $V_0=1, 10, 100$，解方程，程序如下：

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
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
function Poisson2 ()
<p>sys_para = [1    % delta_c<br />
0.5   % gamma<br />
1];  % V</p>
<p>options = bvpset('RelTol',1e-5);%[];    % place holder</p>
<p>x_init = linspace ( 0.0, 25, 5000);</p>
<p>solinit = bvpinit ( x_init, @Poisson2_init);</p>
<p>sol = bvp4c ( @Poisson2_ode, @Poisson2_bc, solinit, options, sys_para);<br />
sol = bvp4c ( @Poisson2_ode, @Poisson2_bc, sol, options, sys_para);</p>
<p>x = sol.x;<br />
y = deval ( sol, x );</p>
<p>close all</p>
<p>plot ( x*0.095,sinh(y(1,:))./(1+2<em>sys_para(2)</em>sinh(y(1,:)/2).^2),'DisplayName','V_0=1', 'Linewidth', 2 );%c_0/pi<br />
legend('-DynamicLegend');<br />
xlabel ( '&lt;--- x/a ---&gt;', 'Fontsize', 16 );<br />
ylabel ( '&lt;--- \rho ---&gt;', 'Fontsize', 16 );<br />
%title ( 'Distribution of positive ion', 'Fontsize', 16 )<br />
grid on<br />
hold all;</p>
<p>sys_para(3) = 10;<br />
sol = bvp4c ( @Poisson2_ode, @Poisson2_bc, sol, options, sys_para);</p>
<p>x = sol.x;<br />
y = deval ( sol, x );</p>
<p>plot ( x*0.095,sinh(y(1,:))./(1+2<em>sys_para(2)</em>sinh(y(1,:)/2).^2),'DisplayName','V_0=10', 'Linewidth', 2);</p>
<p>for i = 20:10:50<br />
sys_para(3) = i;<br />
sol = bvp4c ( @Poisson2_ode, @Poisson2_bc, sol, options, sys_para);<br />
%sol = bvp4c ( @Poisson2_ode, @Poisson2_bc, sol, options, sys_para);<br />
end<br />
for i = 50:5:100<br />
sys_para(3) = i;<br />
sol = bvp4c ( @Poisson2_ode, @Poisson2_bc, sol, options, sys_para);<br />
sol = bvp4c ( @Poisson2_ode, @Poisson2_bc, sol, options, sys_para);<br />
end<br />
%{<br />
sys_para(3) = 100;<br />
sol = bvp4c ( @Poisson2_ode, @Poisson2_bc, sol, options, sys_para);<br />
sol = bvp4c ( @Poisson2_ode, @Poisson2_bc, sol, options, sys_para);<br />
%}</p>
<p>x = sol.x;<br />
y = deval ( sol, x );</p>
<p>plot ( x*0.095,sinh(y(1,:))./(1+2<em>sys_para(2)</em>sinh(y(1,:)/2).^2),'DisplayName','V_0=100', 'Linewidth', 2);<br />
%{<br />
x = sol.x;<br />
y = deval ( sol, x );</p>
<p>close all</p>
<p>plot ( x*0.095,sinh(y(1,:))./(1+2<em>sys_para(2)</em>sinh(y(1,:)/2).^2), 'r-', 'Linewidth', 2 );%c_0/pi<br />
xlabel ( '&lt;--- X ---&gt;', 'Fontsize', 16 );<br />
ylabel ( '&lt;--- Y(X) ---&gt;', 'Fontsize', 16 );<br />
title ( 'Distribution of positive ion', 'Fontsize', 16 )<br />
grid on<br />
filename = 'double-layer';<br />
%figure<br />
%}<br />
return<br />
end</p>
<p>function dydx = Poisson2_ode(x,y,sys_para)</p>
<p>dydx(1) = y(2);<br />
dydx(2) = sinh(y(1))/(1+2<em>sys_para(2)</em>sinh(y(1)/2)^2);</p>
<p>return<br />
end</p>
<p>function bc = Poisson2_bc(ya,yb,sys_para)</p>
<p>bc(1)=ya(1)-sys_para(3);<br />
bc(2)=yb(2);<br />
return<br />
end</p>
<p>function y_init = Poisson2_init(x)</p>
<p>y_init(1)=exp(-x);<br />
y_init(2)=-exp(-x);</p>
<p>return<br />
end</p>
计算结果为

标签: matlab, bvp4c, 离子液体, 泊松-费米方程

添加新评论

captcha
请输入验证码