那看学到什么程度了,微积分说容易也很容易,说困难也很困难.
如果你对微积分的理解就是
\begin{align*} \int_a^b\frac{1}{\sqrt{k^2+x^2}}dx &\underrightarrow{x = ktant}\int_{arctan\frac{k}{a}}^{arctan\frac{k}{b}}\frac{1}{cost}dt\\ &=\int_{arctan\frac{k}{a}}^{arctan\frac{k}{b}}\frac{1}{1-sin^2t}dsint\\ &=\frac{1}{2}ln\frac{1+sint}{1-sint}|_{arctan\frac{k}{a}}^{arctan\frac{k}{b}}\\ &=ln\{\frac{a}{b}\frac{\sqrt{b^2+k^2}+k^2}{\sqrt{a^2+k^2}+k^2}\} \end{align*}
这种简单的积分,那么这种东西一天就可以搞懂.
但是如果你对微积分的理解是
H-H模型
I = C_m\frac{dV_M}{dt} +\overline{g}_Kn^4(V_m-V_K)+\overline{g}_{Na}m^3h(V_m-V_{Na})+\overline{g}_l(V_m-V_l) \\ \frac{dn}{dt} = \alpha_n(V_m)(1-n)-\beta_n(V_m)n \\ \frac{dm}{dt} = \alpha_m(V_m)(1-m) - \beta_m(V_m)m \\ \frac{dh}{dt} = \alpha_h(V_m)(1-h) - \beta_h(V_m)h
这种非线性的体系
在Hodgkin和Huxley的原始论文[1]中，\alpha和\beta函数如下给出\\ \begin{cases} \alpha_n(V_m) = \frac{0.01(10-V)}{exp(\frac{10-V}{10})-1}\\ \alpha_m(V_m) = \frac{0.1(25-V)}{exp(\frac{25-V}{10})-1}\\ \alpha_h(V_m) = 0.07exp(-\frac{V}{20})\\ \beta_n(V_m) = 0.125exp(-\frac{V}{80})\\ \beta_m(V_m) = 4exp(-\frac{V}{18})\\ \beta_h(V_m) = \frac{1}{exp(\frac{30-V}{10})+1}\\ \end{cases}
经典的Hodgkin-Huxley模型将可兴奋细胞的每个部分都当作电路元件来处理（如图所示）。磷脂双分子层表示为电容（C_m）。电压门控离子通道表示为电导（g_n，其中n是特定的离子通道），它依赖于电压和时间。漏通道表示为线性电导（g_L）。驱使离子流动的电化学梯度表示为电压源（E_n），电压源的电压取决于相关离子种类在细胞内和细胞外的浓度的比值。最后，离子泵表示为电流源（i_p）。膜电位表示为V_m。
数学上，流过磷脂双分子层的电流写为 I_c = C_m\frac{dV_m}{dt}
流过给定离子通道的电流为乘积 I_i = g_i(V_m - V_i)
其中 V_i 是第i个离子通道的反转电位。因此，对于具有钠和钾离子通道的细胞，通过细胞膜的总电流为 I = C_m\frac{dV_m}{dt}+g_K(V_m-V_k)+g_{N_a}(V_m-V_{N_a})+g_l(V_m-V_l)
其中I为单位面积的总膜电流， C_m 为单位面积的膜电容， g_K 和 g_{N_a} 分别为单位面积的钾和钠的电导， V_K 和 V_{N_a} 分别为钾和钠的反转电位， g_l 和 V_l 分别为单位面积的漏电导和漏反转电位。这个方程中对时间依赖的项为 V_m 、 g_{N_a} 和 g_K ，其中最后两个电导项也明确地取决于电压。
很明显,我们无法对这四个耦合的微分方程进行解析求解,但是这挡不住我们使用数值解法来求解,这些东西在钱学森的时代已经有了,但是计算机进入普通人的生活还好像21世纪的事情。
我们使用四阶龙格库塔算法实现在持续刺激下电位的变化
我计算了刺激电流 I_{app} 从0到50mV电流的情况
画出如下示意图
各离子在一个周期内的电流情况
可以认为Hodgkin-Huxley模型是一个具有4个状态变量 V_m(t),n(t),m(t),h(t) 的微分方程系统，它们随着时间变化。这个系统很难研究，因为它是一个非线性系统，无法用解析法求解。然而，可以用许多数值方法分析该系统。可以证明某些性质和一般行为（如极限环）是存在的。
我们也可以使用计算机画出来
自然界有很多现象不是由简单的微分方程描述,而是又非线性的方程描述,这就形成了现代物理的一个重要分支,非线性物理,理论物理学家认为非线性是由于量子力学的不确定性原理和经典力学的混沌造成的.我相信这已经不属于钱学森所指的微积分的范畴,但是这些问题又来自于以前微积分无法解决的问题,是一个深奥的领域.
四阶龙格库塔算法实现在持续刺激下电位的变化：
import numpy as np
import matplotlib.pyplot as plt
def alphy_n(u):
return 0.01*(10-u)/(np.exp((10-u)/10)-1)
def alphy_m(u):
return 0.1*(25-u)/(np.exp((25-u)/10)-1)
def alphy_h(u):
return 0.07*np.exp(-u/20)
def beta_n(u):
return 0.125*np.exp(-u/80)
def beta_m(u):
return 4*np.exp(-u/18)
def beta_h(u):
return 1/(np.exp((30-u)/10)+1)
def fV(V,m,h,n):
return (-g_L * (V - E_L) - g_Na * m ** 3 * h * (V - E_Na) - g_K * n ** 4 * (V - E_K) +10)/ C
def fm(V,m,h,n):
return (alphy_m(V) * (1 - m) - beta_m(V) * m)
def fh(V,m,h,n):
return (alphy_h(V) * (1 - h) - beta_h(V) * h)
def fn(V,m,h,n):
return (alphy_n(V) * (1 - n) - beta_n(V) * n)
C = 1
E_Na = 115
g_Na = 120
E_K = -12
g_K = 36
E_L = 10.6
g_L = 0.3
t = np.linspace(0,200,10001)
N = len(t)
dt = t[2] - t[1]
V = np.array([0,0])
m = np.array([0,0.529])
n = np.array([0,0.3177])
h = np.array([0,0.5961])
for i in range(1,N-1):
k1 = fV(V[i],m[i],h[i],n[i])
l1 = fm(V[i],m[i],h[i],n[i])
p1 = fh(V[i],m[i],h[i],n[i])
q1 = fn(V[i],m[i],h[i],n[i])
k2 = fV(V[i]+0.5*dt*k1,m[i]+0.5*dt*l1,h[i]+0.5*dt*p1,n[i]+0.5*dt*q1)
l2 = fm(V[i]+0.5*dt*k1,m[i]+0.5*dt*l1,h[i]+0.5*dt*p1,n[i]+0.5*dt*q1)
p2 = fh(V[i]+0.5*dt*k1,m[i]+0.5*dt*l1,h[i]+0.5*dt*p1,n[i]+0.5*dt*q1)
q2 = fn(V[i]+0.5*dt*k1,m[i]+0.5*dt*l1,h[i]+0.5*dt*p1,n[i]+0.5*dt*q1)
k3 = fV(V[i]+0.5*dt*k2,m[i]+0.5*dt*l2,h[i]+0.5*dt*p2,n[i]+0.5*dt*q2)
l3 = fm(V[i]+0.5*dt*k2,m[i]+0.5*dt*l2,h[i]+0.5*dt*p2,n[i]+0.5*dt*q2)
p3 = fh(V[i]+0.5*dt*k2,m[i]+0.5*dt*l2,h[i]+0.5*dt*p2,n[i]+0.5*dt*q2)
q3 = fn(V[i]+0.5*dt*k2,m[i]+0.5*dt*l2,h[i]+0.5*dt*p2,n[i]+0.5*dt*q2)
k4 = fV(V[i]+dt*k3,m[i]+dt*l3,h[i]+dt*p3,n[i]+dt*q3)
l4 = fm(V[i]+dt*k3,m[i]+dt*l3,h[i]+dt*p3,n[i]+dt*q3)
p4 = fh(V[i]+dt*k3,m[i]+dt*l3,h[i]+dt*p3,n[i]+dt*q3)
q4 = fn(V[i]+dt*k3,m[i]+dt*l3,h[i]+dt*p3,n[i]+dt*q3)
V_next = V[i] + dt/6*(k1+2*k2+2*k3+k4)
m_next = m[i] + dt/6*(l1+2*l2+2*l3+l4)
h_next = h[i] + dt/6*(p1+2*p2+2*p3+p4)
n_next = n[i] + dt/6*(q1+2*q2+2*q3+q4)
V = np.append(V,V_next)
m = np.append(m,m_next)
h = np.append(h,h_next)
n = np.append(n,n_next)
print(m[i],n[i],h[i],V[i])
for i in range(1,N):
V[i] = V[i] - 65
plt.xlabel( 't ms')
plt.ylabel('V mV')
plt.plot(t,V)
plt.show()[1]A quantitative description of membrane current and its application to conduction and excitation in nerve,The Journal of Physiology,Hodgkin AL, Huxley AF,1952
评论区有人说下面的简单,其实是例子简单,我列举一个稍微复杂的.
斑图动力学
在自然界中，斑图现象广泛存在，从动物的皮毛纹理到流体中的对流结构等，这些斑图蕴含着复杂而有序的规律。1952 年，图灵开创性地运用反应 - 扩散系统数学模型描述如虎纹、豹斑等自然斑图，开启了斑图动力学研究的大门。此后，斑图动力学在多个物理领域成为研究热点，其在揭示自然现象背后的物理机制以及推动相关技术发展方面具有关键作用。例如，在流体力学的热对流研究中，通过对瑞利 - 贝纳德对流系统的探索，有助于理解大气环流、海洋热盐对流等宏观现象；在非线性光学系统里，对激光腔光强斑图的研究为光通信、光存储等技术提供了理论支撑；于液晶材料方面，研究液晶分子取向斑图变化对制造高质量显示面板意义重大。本文将深入探讨斑图动力学的模型、解法、图灵斑图特性及线性稳定性分析等核心内容，系统阐述其理论与应用价值。
\[ \begin{cases} \frac{\partial u}{\partial t} = -\frac{1}{\epsilon}u(u-1)[u - (v+b)/a]+\nabla^2u\\ \frac{\partial v}{\partial t} = f(u) -v \end{cases} \]
\[ \begin{cases} f(u) = 0 \text{ if } 0 \leq u < \frac{1}{3}\\ f(u) = 1 - 6.75u(u-1)^2 \text{ if } \frac{1}{3} \leq u \leq 1\\ f(u) = 1 \text{ if } u>1 \end{cases} \]
参数的选择
\[ \begin{cases} a = 0.84\\ b = 0.07\\ \text{如果} 0.01<\epsilon<0.06:\text{可以产生稳定的螺旋波和行波}\\ \text{如果} \epsilon = 0.06 \text{螺旋波将会经历从稳定到不稳定}\\ \text{如果} \epsilon > 0.07 \text{螺旋波破碎,整个系统进入混沌} \end{cases} \]
边界条件
第一种猜测:
零流边界条件: \frac{\partial V}{\partial n} = 0
在非四角边界处:在2D的情况下,用周围4点的值代替中心,但是对于这样的边界周围只有三个点.
如图所示我们可以使用 u_{i,j} 上面的点 u_{i-1,j} ,来代替 u_{i+1,j}
在四角处:
对于 u_{i,j} 我们用 u_{i-1,j} 和 u_{i,j-1} 来代替
但是发现行波并不好
这时我换一种方法:
这时我们发现行波非常完美:
靶波
将中心四点在每一次迭代的时候设置为1,这样可以形成稳定的靶波
螺旋波
要点:产生半列单向传递的行波!在边界处产生一列行波,当行波传递到中间时,将整个上半平面或下半平面设置为0(如果只是设置一列为零无效)
混沌
在模型1中要产生混沌，我们需要令 \epsilon > 0.07
这和产生螺旋波一样，先产生一列行波，然后行波运行到一半时，消除半列波，此时产生的螺旋波会迅速的演化为混沌
消除混沌
当我们在左端加上一个高频的行波,我们可以消除螺旋波.
消除螺旋波
消除螺旋波和混沌的注意事项
如果我们直接消除螺旋波和混沌会导致行波和螺旋波撞不过行波
这时候会被行波撞回去
这时候如果我们选择在产生行波条件发出后调整a和b的值
令a = 0.08,b = 0.06
就可以使行波消除螺旋波和混沌
隐式解法
一维齐次扩散方程矩阵解法的推导
上面我们使用的解决扩散方程的显式解法.
下面我们介绍后向差分法.
对于一维扩散方程 \frac{\partial u(x,t)}{\partial t} = D\frac{\partial^2u(x,t)}{{\partial x}^2}
取 \Delta x = h , \Delta t = \tau 进行离散化,节点坐标为:
\[ \begin{cases} x_i = (i-1)h(i=1,2,\ldots,N)\\ t_k = k\tau(k=1,2,\ldots,k_{max})\\ \end{cases} \]
进过整理 -\xi u_{i+1}^k + (1+2\xi) u_i^k - \xi u_{i-1}^k=u_i^{k-1}
\frac{\partial^2u}{{\partial x}^2} = \frac{u_{i+1}-2u_i+u_{i-1}}{h^2}
可以得到
\frac{u_i^k-u_i^{k-1}}{\tau} = \frac{D}{h^2}(u_{i+1}^k-2u_i^k+u_{i-1}^k)
我们设 \xi = \frac{D \tau}{h^2}
-\xi u_{i+1}^k + (1+2\xi) u_i^k - \xi u_{i-1}^k=u_i^{k-1}
我们设 \xi = \frac{D \tau}{h^2}
-\xi u_{i+1}^k + (1+2\xi) u_i^k - \xi u_{i-1}^k=u_i^{k-1}
\begin{pmatrix} 1 + 2\xi & -\xi & 0 & 0 & 0 \\ -\xi & 1+2\xi & -\xi & 0 & 0 \\ 0 & -\xi & 1+2\xi & 0 & 0 \\ \ldots & \ldots & \ldots & \ldots & \ldots \\ 0 & 0 & 0 & 1+2\xi & -\xi \\ 0 & 0 & 0 &-\xi & 1+2\xi \end{pmatrix} \begin{pmatrix} u_1^k \\ u_2^k \\ u_3^k \\ \ldots \\ u_{n-1}^k \\ u_n^k \end{pmatrix}= \begin{pmatrix} u_1^{k-1} \\ u_2^{k-1} \\ u_3^{k-1} \\ \ldots \\ u_{n-1}^{k-1} \\ u_n^{k-1} \end{pmatrix} \\
一维非齐次的推导
u_t - \nabla^2u = f
\frac{u_i^k-u_i^{k-1}}{\tau} - \frac{u_{i+1}^k-2u_i^k+u_{i-1}^k}{h^2}=f_i^k
u_i^{k-1}=-\frac{\tau}{h^2}u_{i+1}^k+(1+2\frac{\tau}{h^2})u_i^k-\frac{\tau}{h^2}u_{i-1}^k-\tau f_i^k
\begin{pmatrix} f_1^k&1 + 2\xi & -\xi & 0 & 0 \\ f_2^k&-\xi & 1+2\xi & 0 & 0 \\ f_3^k&0 & -\xi & 0 & 0 \\ \ldots&\ldots & \ldots & \ldots & \ldots \\ f_{n-1}^k&0 & 0 & 1+2\xi & -\xi \\ f_n^k&0 & 0 &-\xi & 1+2\xi \end{pmatrix} \begin{pmatrix} 1\\ u_1^k \\ u_2^k \\ u_3^k \\ \ldots \\ u_{n-1}^k \\ u_n^k \end{pmatrix}= \begin{pmatrix} u_1^{k-1} \\ u_2^{k-1} \\ u_3^{k-1} \\ \ldots \\ u_{n-1}^{k-1} \\ u_n^{k-1} \end{pmatrix}
二维非齐次扩散方程的矩阵解法
我们有网格 u_{100\times100}
u_t - D\nabla^2u = f(x,y)
\frac{u_{i,j}^{k+1}-u_{i,j}^k}{\tau} - \frac{D}{h^2}(u_{i-1,j}^k+u_{i+1,j}^k+u_{i,j-1}^k+u_{i,j+1}^k-4u_{i,j})= f_{i,j}
u_{i,j}^{k+1} = u_{i,j}^k+\frac{D\tau}{h^2}(u_{i-1,j}^k+u_{i+1,j}^k+u_{i,j-1}^k+u_{i,j+1}^k-4u_{i,j})+\tau f_{i,j}
设 \frac{D\tau}{h^2} = \xi
u_{i,j}^{k+1} = (1-4\xi)u_{i,j}^k+\xi(u_{i-1,j}^k+u_{i+1,j}^k+u_{i,j-1}^k+u_{i,j+1}^k)+\tau f_{i,j}
对于边界处的特殊处理:
1.四角处的矩阵元(以 u_{11} 为例):
u_{11}^{k+1}=u_{11}^k+\xi(u_{12}^k+u_{21}^k-2u_{11})+\tau f_{i,j}
2.非四角处边界处的矩阵元(以u[:,0]为例)
u_{i,j}^{k+1}+\xi(u_{i-1,j}^k+u_{i+1,j}^k+2u_{i,j+1}^k-4u_{i,j})+\tau f_{i,j} =\xi u_{i,j}^k+\xi(u_{i-1,j}^k+u_{i+1,j}^k+2u_{i,j+1}^k)+\tau f_{i,j}
\begin{pmatrix} \ldots & \ldots & \ldots \\ \ldots & P & \ldots \\ \ldots & \ldots & \ldots \end{pmatrix}_{n^2,n^2+1} \begin{pmatrix} u_{11}^k \\ \ldots \\ u_{n1}^k \\ u_{12}^k \\ \ldots \\ u_{n1}^k \\ \ldots \\ \ldots \\ u_{1n}^k \\ \ldots \\ u_{n,n}^k \end{pmatrix} = \begin{pmatrix} u_{11}^{k+1}\\ \ldots\\ u_{n1}^{k+1}\\ u_{12}^{k+1}\\ \ldots\\ u_{n1}^{k+1}\\ \ldots\\ \ldots\\ u_{1n}^{k+1}\\ \ldots\\ u_{nn}^{k+1}\\ \end{pmatrix}
对于 1<i<n,1<j<n 的情况
u_{i,j}\overrightarrow{}P_{(i-1)n+j,(i-1)n+j} = 1 - 4\xi\\ u_{i,j-1}\overrightarrow{}P_{(i-1)n+j,(i-1)n+j-1} = \xi\\ u_{i,j+1}\overrightarrow{}P_{(i-1)n+j,(i-1)n+j+1} = \xi\\ u_{i-1,j}\overrightarrow{}P_{(i-1)n+j,(i-2)n+j} = \xi\\ u_{i+1,j}\overrightarrow{}P_{(i-1)n+j,in+j} = \xi\\ P_{(i-1)n+j,n^2+1} = f_{ij}\\
P矩阵的前 n^2\times n^2 部分是关于u的系数部分,而最后一列是将 f_{n,n} 排开为 n^2 长度的矢量