4阶泊松-费米方程数值求解

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

\begin{equation} \frac{d^2\phi}{dx^2}-\delta_c^2\frac{d^4\phi}{dx^4}=\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'''(0)=&0 \\ \phi(\infty)=&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
function Poisson4 ()
<p>sys_para = [10    % 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, 50, 1000);</p>
<p>solinit = bvpinit ( x_init, @Poisson4_init);</p>
<p>sol = bvp4c ( @Poisson4_ode, @Poisson4_bc, solinit, options, sys_para);</p>
<p>x = sol.x;<br />
y = deval ( sol, x );<br />
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 ( ' x/a ', 'Fontsize', 16 );<br />
ylabel ( ' \rho ', 'Fontsize', 16 );<br />
title ( 'Dimensionless charge density', 'Fontsize', 16 )<br />
grid on<br />
hold all;</p>
<p>sys_para(3)=10;<br />
sol = bvp4c ( @Poisson4_ode, @Poisson4_bc, sol, options, sys_para);</p>
<p>x = sol.x;<br />
y = deval ( sol, x );<br />
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:5:100<br />
sys_para(3)=i;<br />
sol = bvp4c ( @Poisson4_ode, @Poisson4_bc, sol, options, sys_para);<br />
end</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 />
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 = 'correlation';<br />
%figure<br />
%}<br />
return<br />
end</p>
<p>function dydx = Poisson4_ode(x,y,sys_para)</p>
<p>dydx(1) = y(2);<br />
dydx(2) = y(3);<br />
dydx(3) = y(4);<br />
dydx(4) = (y(3)-sinh(y(1))/(1+2<em>sys_para(2)</em>sinh(y(1)/2)^2))/sys_para(1)^2;</p>
<p>return<br />
end</p>
<p>function bc = Poisson4_bc(ya,yb,sys_para)</p>
<p>bc(1)=ya(1)-sys_para(3);<br />
bc(2)=ya(4);<br />
bc(3)=yb(1);<br />
bc(4)=yb(2);<br />
return<br />
end</p>
<p>function y_init = Poisson4_init(x)</p>
<p>y_init(1)=exp(-x);<br />
y_init(2)=-exp(-x);<br />
y_init(3)=exp(-x);<br />
y_init(4)=-exp(-x);</p>
<p>return<br />
end</p>

计算结果为：

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

添加新评论

captcha
请输入验证码