计算科学 - 各向同性向列相变的计算 - 吾爱随笔录

各向同性向列相变的计算

计算科学 C 牛顿法

2021-12-29 20:13:18

本文讨论了各向同性-向列相变背后的理论。此外，给出了一种算法来计算这种相变的一些性质。

我用 C 语言编写了一个程序，它为我提供了给定浓度的自由能、压力、化学势和向列序参数 $c$ 系统中的粒子。现在我想找到浓度 $c_I$ 和 $c_N$ ，它定义了系统相分离的间隔，即它部分是向列的，部分是各向同性的，以降低其自由能。

正如文章所解释的，需要求解方程组

p_{iso} (c_{I}) = p (c_{N}) μ_{iso} (c_{I}) = μ (c_{N}),

$p_{\text{iso}}(c_I) = p(c_N) \\ \mu_{\text{iso}}(c_I) = \mu(c_N),$ 为了

c_{I, N}

$c_{I,N}$ . 表达式

p_{iso}

$p_{\text{iso}}$ 和

μ_{iso}

$\mu_{\text{iso}}$ 是分析已知的，我可以计算

p

$p$ 和

μ

$\mu$ 对于浓度的某些离散值。

我的问题是：如何在 C 中求解这两个方程？我知道可以使用 Newton-Rhapson 以数值方式求解函数的零点。但是，因为我只知道功能 $p$ 和 $\mu$ 对于某些值 $c$ ，我真的不知道如何解决这个问题。

有人可以帮我找到一种方法吗？（最好是某种算法。）免责声明：我开始学习 C，但我不是很擅长。

1个回答

我没有足够的声誉来写评论，所以这个作为答案。

我认为您已经编写了解决角度分布的代码 $\psi$ 给定浓度的数值 $c$ ，使用论文中概述的自洽程序。现在您可以将自由能写为 phi 的函数，而压力是 $c$ -自由能的导数，原则上，你可以得到任何值的压力 $c$ .

我有点困惑，在推导压力时，似乎 $\psi$ s 不依赖于 $c$ . 他们显然做到了。所以我相信这里有一个隐含的近似值。如果保持相同的近似值，则可以手动写下压力的导数。然后，您可以在任何给定的情况下评估压力及其导数 $c$ （这就是为什么当你说你只知道它们的某些值时我不确定你的意思 $c$ ）。现在你有了这些函数（化学势也是如此），写下牛顿-拉夫森方案应该没有问题（你有一个两个变量，两个函数设置，所以你必须区分两个函数相对于每个变量，然后反转得到的雅可比矩阵）。

编辑： 所以我决定自己用“Python”编写代码。我说“Python”而不是 Python，因为我将它用作工具箱中的一个快速下拉菜单。我将 ipython 与 -pylab 开关一起使用，因此它将加载部分 scipy、numpy 进行计算。

Ntheta  = 128
ks = arange(Ntheta+1)[1:]

thetas = pi/2.*ks/(Ntheta+1)

deltas = zeros_like(thetas)
deltas[0]    = 1 - (cos(thetas[0]) + cos(thetas[1]))/2.
deltas[1:-1] = (cos(thetas[:-2]) - cos(thetas[2:]))/2.
deltas[-1]   = (cos(thetas[-1]) + cos(thetas[-2]))/2.

Nphi = 256
js = arange(Nphi+1)[1:]
phis = 2.*pi*js/(Nphi+1)
Ks = zeros((Ntheta, Ntheta))
for k in range(Ntheta):
  for l in range(Ntheta):
    g = sqrt(1. - (cos(thetas[k])*cos(thetas[l]) + sin(thetas[k])*sin(thetas[l])*cos(phis))**2)
    Ks[k, l] = 2.*pi/(Nphi+1) * (3./2*g[0] + sum(g[1:-1]) + 3./2*g[-1])

def thefunc(c):
  taus = (c/pi)**2 * exp(-2*c**2*thetas**2/pi)

  for i in range(128):
    As = 16.*c/pi * dot(Ks, deltas*taus)
    Z  = 4.*pi*sum(deltas*exp(-As))
    psis = 1./Z * exp(-As)
    taus = psis.copy()

  S = 4.*pi*sum(deltas*psis*(3.*cos(thetas)**2-1.)/2.)
  theint = dot(deltas*psis, dot(Ks, deltas*psis))
  p = c + 32.*c**2*theint
  pder = 1. + 64.*c*theint
  f  = log(c) - 1. + 4.*pi*sum(deltas*psis*log(psis)) + 32.*c*theint
  mu = f + p/c
  muder = pder/c
  return p, pder, mu, muder

def isofunc(c):
  p = c + c**2
  pder = 1. + 2.*c
  mu = log(c/(4.*pi)) + 2.*c
  muder = 1./c + 2.
  return p, pder, mu, muder

def totalfunc(ct, c):
   pt, pdert, mut, mudert = thefunc(ct)
   p, pder, mu, muder     = isofunc(c)
   F = -array([pt - p, mut - mu])
   J = array([[pdert, -pder], [mudert, -muder]])
   dc = solve(J, F)
   return ct+dc[0], c+dc[1]

ct = 4.
c  = 3.
for i in range(100):
  ct, c = totalfunc(ct, c)

这里thefunc计算自洽解 $\psi$ 返回压力，化学势及其衍生物。isofunc对各向同性相也是如此。totalfunc进行 Newton-Raphson 迭代。很抱歉函数的名称。不过，这些变量应该是不言自明的。

其它你可能感兴趣的问题

上一篇Matlab 表示图下一篇图像边界上折叠矩阵的行为（高斯模糊）