计算物理课程笔记
第一章基本数值运算
微分
向前 \[ \begin{aligned} &f(x_0+h)=f(x_0)+hf'(x_0)+\frac{h^2}{2}f''(\xi)\\ &\frac{f(x_0+h)-f(x_0)}{h}=f'(x_0)+\frac{h}{2}f''(\xi)\\ &R=-\frac{h}{2}f''(\xi) \end{aligned} \]
向后 \[ \begin{aligned} &f(x_0-h)=f(x_0)-hf'(x_0)+\frac{h^2}{2}f''(\xi)\\ &\frac{f(x_0)-f(x_0-h)}{h}=f'(x_0)-\frac{h}{2}f''(\xi)\\ &R=\frac{h}{2}f''(\xi) \end{aligned} \]
中心差分 \[ \begin{aligned} &f(x_0+h)=f(x_0)+hf'(x_0)+\frac{h^2}{2}f''(x)+\frac{h^3}{3!}f'''(\xi_{+})\\ &f(x_0-h)=f(x_0)-hf'(x_0)+\frac{h^2}{2}f''(x)-\frac{h^3}{3!}f'''(\xi_{-})\\ &\frac{f(x_0+h)-f(x_0-h)}{2h}=f'(x_0)+\frac{h^2}{12}(f'''(\xi_{+})+f'''(\xi_{-}))\\ &R =-\frac{h^2}{12}(f'''(\xi_{+})+f'''(\xi_{-})) \end{aligned} \]
五点公式 \[ \begin{aligned} &f(x_0+h)=f(x_0)+hf'(x_0)+\frac{h^2}{2}f''(x)+\frac{h^3}{3!}f'''(x)\\ &\quad\quad\quad\quad\quad +\frac{h^4}{4!}f^{(4)}(x)+\frac{h^5}{5!}f^{(5)}(\xi_{+})\\ & f(x_0-h)=f(x_0)-hf'(x_0)+\frac{h^2}{2}f''(x)-\frac{h^3}{3!}f'''(x)\\ &\quad\quad\quad\quad\quad +\frac{h^4}{4!}f^{(4)}(x)-\frac{h^5}{5!}f^{(5)}(\xi_{-})\\ &f(x_0+2h)=f(x_0)+2hf'(x_0)+2h^2f''(x)+\frac{4h^3}{3}f'''(x)\\ &\quad\quad\quad\quad\quad+\frac{2h^4}{3}f^{(4)}(x)+\frac{(2h)^5}{5!}f^{(5)}(\xi_{+})\\ &f(x_0-2h)=f(x_0)-2hf'(x_0)+2h^2f''(x)-\frac{4h^3}{3}f'''(x)\\ &\quad\quad\quad\quad\quad+\frac{2h^4}{3}f^{(4)}(x)-\frac{(2h)^5}{5!}f^{(5)}(\xi_{-})\\ & \frac{8(f(x_0+h)-f(x_0-h))-(f(x_0+2h)-f(x_0-2h))}{12h}\\ &\quad\quad\quad\quad\quad=f'(x)+O(h^5)\\ & R = O(h^5) \end{aligned} \]
积分
langrange插值法
一阶插值 函数已知(\(x_0\),\(y_0\)),(\(x_1\),\(y_1\))两点。可以取函数。 \[ y = y_0\frac{x-x_1}{x_0-x_1}+y_1\frac{x-x_0}{x_1-x_0} \] 实际上可以引入lagrange插值基函数 \[ l_0(x)=\frac{x-x_1}{x_0-x_1}\\ l_1(x)=\frac{x-x_0}{x_1-x_0} \] 插值基函数的特点:再角标对应的点数值为1，在其他角标对应的点数值为0
高阶插值(比如说n阶)
函数已知点(\(x_0\),\(y_0\)),(\(x_1\),\(y_1\)) ... (\(x_n\),\(y_n\))
可以取函数: \[ y = y_0l_0(x)+y_1l_1(x)+y_2l_2(x)...+y_nl_n(x) \] 插值基函数: \[ l_i(x)=\frac{(x-x_0)(x-x_1)...(x-x_{i-1})(x-x_{i+1})...(x-x_n)}{(x_i-x_0)(x_i-x_1)...(x_i-x_{i-1})(x_i-x_{i+1})...(x_i-x_n)} \] 比如说二阶插值 \[ y = y_0\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}+y_1\frac{(x-x_0)(x-x_1)}{(x_1-x_0)(x_1-x_2)}+y_2\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)} \] 比如说三阶插值 \[ \begin{aligned} y = y_0\frac{(x-x_1)(x-x_2)(x-x_3)}{(x_0-x_1)(x_0-x_2)(x_0-x_3)}+y_1\frac{(x-x_0)(x-x_2)(x-x_3)}{(x_1-x_0)(x_1-x_2)(x_1-x_3)}+\\y_2\frac{(x-x_0)(x-x_1)(x-x_3)}{(x_2-x_0)(x_2-x_1)(x_2-x_3)}+y_3\frac{(x-x_0)(x-x_1)(x-x_2)}{(x_3-x_0)(x_3-x_1)(x_3-x_2)} \end{aligned} \]
runge现象 插值阶数比较高时反而效果不好
三种插值积分方式
辛普生成公式(simpson)
考虑三个点(\(x_{i}\),\(y_{i}\)),(\(x_{i+1}\),\(y_{i+1}\)),(\(x_{i+2}\),\(y_{i+2}\)) \[ y = y_{i}\frac{(x-x_{i+1})(x-x_{i+2})}{(x_i-x_{i+1})(x_i-x_{i+2})}+y_{i+1}\frac{(x-x_i)(x-x_{i+2})}{(x_{i+1}-x_{i})(x_{i+1}-x_{i+2})}\\ +y_{i+2}\frac{(x-x_i)(x-x_{i+1})}{(x_{i+2}-x_i)(x_{i+2}-x_{i+1})} \] 将y从\(x_i\)积分到\(x_{i+2}\) \[ I = \int_{x_{i}}^{x_{i+2}}ydx \] 积分得到(草稿E1) \[ \begin{aligned} &I = y_i\frac{-\frac{1}{3}(x_i-x_{i+2})^2+x_i(x_i-x_{i+2})+\frac{1}{2}(x_{i+1}+x_{i+2})(x_i+x_{i+2})-x_i^2-x_{i+1}x_{i+2}}{x_i-x_{i+1}}\\ & +y_{i+1}\frac{1}{(x_{i+1}-x_{i})(x_{i+1}-x_{i+2})}\Big(\frac{1}{3}\left((x_{i+2}-x_{i+1}\right)^3-(x_i-x_{i+1})^3) \\ &\quad\quad +\frac{1}{2}\left((x_{i+2}-x_{i+1})^2-(x_i-x_{i+1})^2\right)(2x_{i+1}-x_i-x_{i+2})\\ &\quad\quad+(x_{i+1}-x_i)(x_{i+1}-x_{i+2})(x_{i+2}-x_i)\Big)\\ &+y_{i+2}\frac{1}{(x_{i+2}-x_i)(x_{i+2}-x_{i+1})}\Big(\frac{1}{3}\left((x_{i+2}-x_{i+1})^3-(x_i-x_{i+1})^3\right)\\ &\quad\quad +\frac{1}{2}(x_{i+1}-x_i)\left((x_{i+2}-x_{i+1})^2-(x_i-x_{i+1})^2\right)\Big) \end{aligned} \] 若认为\(x_i\),\(x_{i+1}\),\(x_{i+2}\)之间相差h。则简化积分为 \[ I=\frac{1}{3}h(y_i+4y_{i+1}+y_{i+2}) \] 在此情况下，若将积分区间划分为n等分(从\(x_0\)到\(x_n\),n一定是2的整数倍),h表示相邻两点之间的距离。可以将积分表示为 \[ I = \frac{1}{3}h(f(x_0)+4\Sigma_{j=0}^{\frac{n}{2}-1}f(x_0+h+2jh)+2\Sigma_{i=0}^{\frac{n}{2}-2}f(x_0+2h+2jh)+f(x_n)) \]
simpson\(\frac{3}{8}\)算法。Newton积分公式
如果已知四个点(\(x_{i}\),\(y_{i}\)),(\(x_{i+1}\),\(y_{i+1}\)),(\(x_{i+2}\),\(y_{i+2}\)),(\(x_{i+3}\),\(y_{i+3}\)) \[ \begin{aligned} y =& y_i\frac{(x-x_{i+1})(x-x_{i+2})(x-x_{i+3})}{(x_i-x_{i+1})(x_i-x_{i+2})(x_i-x_{i+3})}\\ &+y_{i+1}\frac{(x-x_i)(x-x_{i+2})(x-x_{i+3})}{(x_{i+1}-x_i)(x_{i+1}-x_{i+2})(x_{i+1}-x_{i+3})}\\ &+y_{i+2}\frac{(x-x_i)(x-x_{i+1})(x-x_{i+3})}{(x_{i+2}-x_i)(x_{i+2}-x_{i+1})(x_{i+2}-x_{i+3})}\\ &+y_{i+3}\frac{(x-x_i)(x-x_{i+1})(x-x_{i+2})}{(x_{i+3}-x_i)(x_{i+3}-x_{i+1})(x_{i+3}-x_{i+2})} \end{aligned} \] 在区间\(x_{i}\)到\(x_{i+3}\)之间积分得(认为间距相等且是h) \[ I =\frac{3}{8}h(f(x_i)+3f(x_{i+1})+3f(x_{i+2})+f(x_{i+3})) \] 若将积分区间分为从\(x_0\)到\(x_n\)一共n个区间,且n是3得整数倍。 \[ \begin{aligned} I = \frac{3}{8}h(f(x_0)+\Sigma_{j=0}^{j=\frac{n}{3}-1}3(f(x_0+h+3jh)+f(x_0+2h+3jh))\\+ \Sigma_{j=0}^{j=\frac{n}{3}-2}2f(x_0+3h+3jh)+f(x_0+nh)) \end{aligned} \]
线性插值积分
如果已知两个点的坐标(\(x_i\),\(y_i\)),(\(x_{i+1}\),\(y_{i+1}\)) \[ I = \frac{1}{2}h(f(x_i)+f(x_{i+1})) \] 若将区间划分为n等分，从\(x_0\)到\(x_n\)。每一份之间距离为h。 \[ I = \frac{1}{2}h(f(x_0)+\Sigma_{j=0}^{j=n-2}2f(x_0+h+jh)+f(x_0+nh)) \]
有一个二重积分的例子几乎用到了所有的技巧:
1 | "F:\BaiduSyncdisk\计算物理与程序设计\第二次课PPT及源程序\第二次课\chap1_integration_double.m" |
寻找函数的根
二分法
二分法的流程:初始化f(x),a,b,\(\epsilon\)(误差标准)。取x = (a+b)/2。若|b-a|<\(\epsilon\),则结束运算,输出x。计算f(x),若f(a)f(x)>0,则a = x,负责令b = x,返回x=(a+b)/2。
一个老师的例子(\(f(x)=2x^3-5x-1\)),matlab 技巧方面,用了二分法(封装在函数里面)。
1 | "F:\BaiduSyncdisk\计算物理与程序设计\第三次课PPT及源程序\chap1_example_3_bisection_function.m" |
newton-Raphson方法
newton-Raphson方法的流程: 初始化\(x_0\),误差标准\(\delta\),\(\epsilon\),设置迭代次数k=0。计算\(x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}\),迭代k = k+1。若\(|f(x_k)|\leq\delta\land|x_{k+1}-x_k|<\epsilon\)则退出循环。
老师的一个例子用牛顿法求解\(f(x)=e^x-1.5-tan^{-1}x\)初始点\(x_0=-7.0\),又用二分法在[-16,-7]上面的解
1 | "F:\BaiduSyncdisk\计算物理与程序设计\第三次课PPT及源程序\chap1_example_4_bisection_newton.m" |
弦割法
需要两个启动点 \[ f'(x_k) = \frac{f(x_k)-f(x_{k-1})}{x_k-x_{k-1}}\\ x_{k+1} = x_{k}-\frac{f(x_k)}{f'(x_k)}\\ \]
停止的标志: \[ |x_{k+1}-x_{k}|<\epsilon \land |f(x_k)|<\delta \] 弦割法的计算流程。初始化\(x_0\),\(x_1\),误差标准\(\delta\)和\(\epsilon\),设置k=0。若\(|f(x_K)|\leq\delta\land|x_K-x_{k+1}|\leq\epsilon\)则停止。计算\(x_{k+1}=x_k-f(x_K)\frac{x_k-x_{k-1}}{f(x_k)-f(x_{k-1})}\)。k=k+1,转到2。
老师的例:\(xe^x-1=0,x_0=0.5,x_1=0.6\),例子程序。
1 | "F:\BaiduSyncdisk\计算物理与程序设计\第三次课PPT及源程序\chap1_example_5_bisection_newton_secant.m" |
- 三种求根方法的总结 最好先用二分法找到大概位置，再用另外两个方法找到根的具体位置。
方程多根
先逐步搜索找到根存在区间，再用上面三种方法找到精确值。 老师的一个例子(程序比较复杂233)
1
"F:\BaiduSyncdisk\计算物理与程序设计\第三次课PPT及源程序\chap1_example_6_several.m"
Matlab程序
- fzero(f,[3,4]) (f是句柄函数)
- roots(c) 多项式的全部零点
- fsolve 非线性方程组的数值解
微分方程求解
初值问题
- 初值问题的题目形式 \[ \frac{dy}{dx}=f(x,y)\\ x\in(a,b)\\ y(x=a)=y_0 \]
Euler 法
\[ y_{n+1}=y_{n}+hf(x_n,y_n) \]
- 老师课上的例子 \[
m\frac{d^2x}{d^2t}=f(x,t)\\
p = m\frac{dx}{dt}\\
\frac{dp}{dt}=f(x,t)\\
\frac{dx}{dt}=p/m\\
x_{n+1}=hp_{n}/m+x_{n+1}\\
p_{n+1}=hf(x_n,tn)+p_n
\]
1
"F:\BaiduSyncdisk\计算物理与程序设计\第四次课PPT及源程序\第四次课PPT及源程序\chapter2_example_2_harmonic_oscillator.m"
Taylor 级数法
误差大概是\(h^2\)的量级。 \[ y(x_n+h)=y(x_n)+hy'(x_n)+\frac{1}{2}h^2y''(x_n) \] 实际上 \[ y'(x_n)=f(x_n,y_n)\\ y''(x_n)=\frac{\partial }{\partial x}f(x_n,y_n)+f\frac{\partial}{\partial y}f(x_n,y_n) \]
Adams-Bashforth 二步法
需要多个启动点的都叫多步法。 \[ y_{n+1}=y_n+\int_{x_n}^{x_{n+1}}f(x,y(x))dx\\ y_{n+1} =y_n+h(\frac{3}{2}f_n-\frac{1}{2}f_{n-1}) \]
AB四步法
\[ y_{n+1}=y_n+h\frac{1}{24}(55f_n-59f_{n-1}+37f_{n-2}-9f_{n-3}) \]
启动点用Taylor级数法得到。
- 上课的一个例题 \[
\frac{dy}{dx}=-xy y(0)=1
\]
1
"F:\BaiduSyncdisk\计算物理与程序设计\第五次课-计算物理第2章-常 微分方程的初值问题\chapter2_example_4_Adams_Bashforth.m"
隐式法
类似于AB方法，不过在插值时用了\(f_{n+1}\)这一项。例如:
因为方程左右都含有\(y_{n+1}\),所以叫做隐式法。
Adams-Moulton一步法 \[ \begin{aligned} y(x_{n+1})&=y(x_n)+\int_{x_n}^{x_{n+1}}f(x,y(x))dx\\ &=y(x_n)+\int_{x_n}^{x_{n+1}}\left(f_n\frac{x-x_{n+1}}{x_n-x_{n+1}}+f_{n+1}\frac{x-x_n}{x_{n+1}-x_n}\right)dx\\ &=y_n+\frac{h}{2}[f_n+f_{n+1}]+O(h^3) \end{aligned} \]
Adams-Moulton二步法 \[ y_{n+1}=y_n+\frac{h}{12}(5f_{n+1}+8f_n-f_{n-1}) \]
AM三步法 \[ y_{n+1}=y_n+\frac{h}{24}(9f_{n+1}+19f_n-5f_{n-1}+f_{n-2}) \]
对隐式法,AB法,Euler法的一点点讨论
AM隐式法不能直接使用,但是AB和Euler 方法可以直接使用,尔后用AM法进行矫正。
有一个常用的这样的算法是AB四步法和AM三步法的结合
改进Euler方法-梯形公式
\[ y_{n+1}=y_n+hf_n\\ y_{n+1}=y_n+h\frac{1}{2}(f_{n+1}+f_n) \]
第一步是预估,第二步是矫正
Rung-kutta算法
基本思想 \[ y_{n+1}=y_n+h\Sigma_{i=1}^{N}\lambda_iK_1\\ K_1=f(x_n,y_n)\\ K_i=f(x_n+hc_i,y_n+hc_i(\Sigma_{j=1}^{i-1}a_{ij}K_j)) \] 对上面的式子进行taylor展开,要求尽可能多的项数和级数相同。
二阶 \[ y_{n+1}=y_n+\frac{h}{2}K_1+\frac{h}{2}K_2\\ K_1=f(x_n,y_n)\\ K_2=f(x_n+h,y_n+hK_1) \]
三阶 \[ \begin{aligned} &y_{n+1}=y_n+\frac{h}{6}(K_1+4K_2+K_3)\\ &K_1=f(x_n,y_n)\\ &K_2=f(x_n+\frac{h}{2},y_n+h\frac{K_1}{2})\\ &K_3=f(x_n+h,y_n+h(-K_1+2K_2)) \end{aligned} \]
四阶(认为是最常用的方法) \[ \begin{aligned} &y_{n+1}=y_n+\frac{h}{6}(K_1+2K_2+2K_3+K_4)\\ &K_1=f(x_n,y_n)\\ &K_2=f(x_n+\frac{h}{2},y_n+\frac{h}{2}K_1)\\ &K_3=f(x_n+\frac{h}{2},y_n+\frac{h}{2}K_2)\\ &K_4=f(x_n+h,y_n+hK_3) \end{aligned} \]
注:对于多因变量(比如\(y_1,y_2,y_3\)),是类似的。
稳定性问题
取实验方程 \[ y'=\lambda y \] 假设在\(y_n\)引入了误差\(\rho_n\),如果在后面误差绝对值不增加,就说这个数值方法对步长h和\(\lambda\)常数是稳定的。步长h和\(\lambda\)有其相应的允许范围,叫做该方法的绝对稳定区域。
- 对于四阶RK算法\(\frac{dy}{dx}=\lambda y\) \[ \rho_{n+1}=[1+(h\lambda)+\frac{(h\lambda)^2}{2!}+\frac{(h\lambda)^3}{3!}+\frac{(h\lambda)^4}{4!}]\rho_n \] 容易得到稳定区域是\(\lambda h\in(-2.78,0)\) ###### 常微分方程的numerov算法
这个算法只能用于特定的微分方程。(很多时候K写为了\(k^2\))
\[
\frac{d^2y}{dx^2}+K(x)y=S(x)
\] 利用二阶导数的差分公式 \[
\frac{y_{n+1}+y_{n-1}-2y_n}{h^2}=\frac{2\frac{1}{2}y^{(2)}h^2+O(h^4)}{h^2}=\frac{2\frac{1}{2}y^{(2)}h^2+2\frac{1}{4!}y^{(4)}h^4+O(h^6)}{h^2}
\] 对微分方程求导: \[
\begin{aligned}
\frac{d^4y}{dx^4}&=\frac{d^2(S(x)-K(x)y)}{dx^2}\\&
=\frac{S_{n+1}+S_{n-1}-2S_n}{h^2}-\frac{K_{n+1}y_{n+1}+K_{n-1}y_{n-1}-2K_{n}y_n}{h^2}+O(h^4)
\end{aligned}
\] 而: \[
\begin{aligned}
S_n-K_ny_n=&\frac{y_{n+1}+y_{n-1}-2y_n}{h^2}\\
&-\frac{1}{12}(S_{n+1}+S_{n-1}-2S_n-K_{n+1}y_{n+1}-K_{n-1}y_{n-1}+2K_{n}y_n)\\
&+O(h^4)
\end{aligned}
\] 可以将上式写为\(y_{n+1},y_{n-1},y_n\)的递推关系 \[
\begin{aligned}
y_{n+1}&(1+\frac{1}{12}h^2K_{n+1})+y_{n-1}(1+\frac{1}{12}h^2K_{n-1})+2y_n(-1+\frac{5}{12}h^2K_n)\\
&=\frac{h^2}{12}(S_{n+1}+10S_n+S_{n-1})
\end{aligned}
\] 这个解法不可以自启动,一般设置启动点的方法是泰勒级数法。 \[
y_1=y_0+hy_0'+\frac{1}{2}h^2y_0''
\]
第三章-常微分方程的边值问题与本征值问题
边值问题的解
线性边值问题的迭加法
问题描述 \[ y''=f(x,y,y') \] 如果可以写为 \[ y''+p(x)y'+q(x)y=f(x) \] 则方程组是线性的。
如果边值问题可以被描述为\(y(x=a)=\alpha,y(x=b)=\beta\),则可以视为两组边值结果线性叠加的效果\((y(a)=\alpha,y'(a)=0,y''+p(x)y'+q(x)y=f(x));(y(a)=0,y'(a)=1,y''+p(x)y'+q(x)y=0)\)
将他们的解分别记为\(y_1\)与\(y_2\),则总的解可以写为。 \[ y(x)=y_1(x)+\frac{\beta-y_1(b)}{y_2(b)}y_2(x) \]
非线性边值问题的打靶法
定义初始条件:\(y(a)=\alpha,y'(a)=s\),则可以根据常微分方程的初值问题的解法得到\(y(b,s)\)
考虑方程 \[ F(s)=y(b,s)-\beta \] 求下一个点的公式可以用弦割法: \[ s_{k+1} = s_{k}-\frac{s_k-s_{k-1}}{y(b,s_k)-y(b,s_{k-1})}(y(b,s)-\beta) \]
- 给定初始点的斜率猜测值S1；
- 用常微分方程的初值问题解法（如RK算法）求解y(b,S1);
- 给出另一个斜率猜测值S2作为弦割法的第二个启动点;
- 用常微分方程的初值问题解法（如RK算法）求解y(b,S2);
- IF (abs(y(b,S1)-y(b))<精度 ), y(x)=y(x,S1);
- IF (abs(y(b,S2)-y(b))<精度 ), y(x)=y(x,S2); else
- 利用弦割法迭代公式求出下一个斜率S3；
- 用常微分方程的初值问题解法（如RK算法）求解y(b,S3);
- IF (abs(y(b,S3)-y(b))<精度 ), y(x)=y(x,S3); else S1=S2; S2=S3；回到第7步
- 得到边值问题的解y(x)=y(x,S3);
- 三类边值问题对应打靶法的参数选取
- 第一类边值问题\(y(a)=\alpha,y(b)=\beta\),显而易见,取\(y'(a)=s\)为参数
- 第二类边值问题\(y'(a)=\alpha,y'(b)=\beta\),取\(y(a)\)为参数
- 第三类边值问题\(y'(a)=\beta_0+\alpha_0y(a),y'(b)=\beta_1+\alpha_1y(b)\)可以将\(y(a)\)取为参数\(y_0\),则\(y'(a)=\beta_0+\alpha_0y_0\)
线性常微分方程的差分方法
\[y''+p(x)y'+q(x)y=f(x)\]
将差分表示的微商带入 \[ \frac{y_{n+1}+y_{n-1}-2y_n}{h^2}+p(x_n)\frac{y_{n+1}-y_{n-1}}{2h}+q(x_n)y_n=f(x) \] 即 \[ y_{n-1}(2-hp(x_n))+y_n(-4+2h^2q(x_n))+y_{n+1}(2+hp(x_n))=2h^2f(x_n) \]
考虑 \(y_0,y_1...y_{n-1}\) (n-1)个点,与n-1个方程
\[ \begin{equation} \begin{aligned} \begin{bmatrix} 2-hp_1 & -4+2h^2q_1 & 2+hp_1 & ... & ... & ... \\ ... & 2-hp_2 & -4+2h^2q_2 & 2+hp_2 & ... & ... \\ ... & ... & ... & ... & ... & ... \\ ... & ... & ... & ... & ... & ... \\ ... & ... &... & 2-hp_{n-1} & -4+2h^2q_{n-1} & 2-hp_{n-1} \\ \end{bmatrix} \begin{bmatrix} y_0\\ y_1\\ ...\\ y_{n-1}\\ y_{n} \end{bmatrix}\\ = \begin{bmatrix} 2h^2f_0\\ 2h^2f_1\\ ...\\ 2h^2f_{n-1}\\ 2h^2f_{n} \end{bmatrix} \end{aligned} \end{equation} \]
但是有n+1个未知数,需要再带上两个边界的方程。总之就是把上面的矩阵扩展为方形。
本征值问题的解:
比如说波动方程的分离变量后的空间部分 \[ y''-k^2y=0 \] 可以取初始条件\(y(x=0)=\alpha,y'(x=0)=\delta\).\(\delta\)实际上可以随便取值.....因为方程是其次的,而且在边界点的函数值刚好是0。然后就用二分法的打靶法,将k设置为参量,即可。 ##### 求解薛定谔方程的本征值问题 \(schr\ddot{o}dinger\) -Equation一般表示为: \[ \begin{align*} \frac{\partial^2}{\partial x^2}\psi(x)+\frac{2m}{\hbar^2}(E-V(x))\psi=0\\ \frac{\partial^2}{\partial x^2}\psi(x)+(\epsilon-v(x))\psi=0 \end{align*} \] 其中,\(\epsilon = \frac{2m}{\hbar^2}E\), \(v(x)=\frac{2m}{\hbar^2}V(x)\). 一般认为这类方程只在一定的区间内需要被考虑,比如说\([x_{min},x_{max}]\),在这个区间之外,一般认为\(v(x)\to +\infty\) , 所以相当于波函数就是0了。而在这个区间内,我们经常把势能写成一个凹函数,并且满足\(v(x_{min})=v(x_{max})=0\)。我们一般考虑束缚态,这样E<0。这样，满足\(E-V(x)\)的区间我们可以标记为\([x_{m1},x_{m2}]\)。 我们的解法是分别给定初值条件: \[ \begin{align} \begin{cases} \psi_{<}(x_{min})=0\\ \frac{\partial}{\partial x}\psi_{<}(x_{min})=\delta \end{cases} \\ \begin{cases} \psi_{>}(x_{\max})=0\\ \frac{\partial}{\partial x}\psi_{>}(x_{max})=\delta' \end{cases} \end{align} \] 这样可以得到两个解\(\psi_{<}(x)\)和\(\psi_{>}(x)\)。我们要做的具体工作是: 先设定本征值,解出两个解,再对其中一个乘以系数,使得他们在\(\psi(x_{m1})\)处的函数值相等,然后再用向后差分的方法,对比两者再\(x_{m1}\)点的一阶导数是否一样,如果一样,就认为这个本征值就是方程的本征值。
第四章 偏微分方程的数值求解
偏微分方程的分类
二阶偏微分方程常见的形式可以写为 \[ AU_{xx}+2BU_{xy}+CU_{yy}+DU_x+EU_y+FU=0 \] 其中,\(A,B,C\)为参数且取决于x,y。 如果再xy平面上有\(A^2+B^2+C^2>0\) ,则该偏微分方程再该平面上为二阶偏微分方程。 * 如果\(B^2-AC<0\), 则为椭圆形方程。 * 如果\(B^2-AC=0\), 则为抛物型方程。 * 如果\(B^2-AC>0\), 则为双曲型方程。 #### 椭圆形方程的求解 ##### 什么是椭圆形方程 Poisson方程在这里实际上是椭圆形方程的代表 Possion方程可以写为: \[ \Delta u = \nabla^2(u)=\frac{\partial^2u}{\partial x^2} +\frac{\partial^2u}{\partial y^2}=S(x,y) \] ##### Jacobi 迭代法 将空间划分为一些格点\(u_{i,j}\)用来表示x方向第i个, y方向第j个格点。 先假设u有一个分布叫做\(u^0\), 利用差分法,迭代得到下一次的分布。 \[ \frac{u_{i+1,j}^n+u_{i-1,j}^n-2u_{i,j}^{n+1}}{h_1^2}+\frac{u_{i,j+1}^n+u_{i,j-1}^n-2u_{i,j}^{n+1}}{h_2^2}=S_{i,j}^n \] 上面的式子中,n代表上一次的分布,n+1代表下一次的分布,理论上,当迭代前后在一个格点上函数值没有变化时，说明找到了possion方程的解。 将上面的式子化简之后得到的式子是: \[ u_{i,j}^{n+1}=(\frac{u_{i+1,j}^n+u_{i-1,j}^n}{h_1^2}+\frac{u_{i,j+1}^n+u_{i,j-1}^n}{h_2^2}-S_{i,j}^n)\frac{1}{\frac{2}{h_1^2}+\frac{2}{h_2^2}}\tag{1} \] ##### Gauss-Seidel松弛迭代法 在[[计算物理#Jacobi 迭代法|Jacobi迭代法]],我们在迭代生成\(u^{n+1}_{i,j}\)时,实际上已经有算出\(u_{i-1,j}^{n+1},u_{i,j-1}^{n+1}\)。在这样的情况下,我们可以用已经算出来的函数值来修正Jacobi迭代法。 按照这个思想, 迭代后的函数值可以表示为: \[ \begin{aligned} u_{i,j}^{n+1}=(1-w)(\frac{u_{i+1,j}^n+u_{i-1,j}^n}{h_1^2}+\frac{u_{i,j+1}^n+u_{i,j-1}^n}{h_2^2}-S_{i,j}^n)\frac{1}{\frac{2}{h_1^2}+\frac{2}{h_2^2}}+\\ w(\frac{u_{i+1,j}^n+u_{i-1,j}^{n+1}}{h_1^2}+\frac{u_{i,j+1}^n+u_{i,j-1}^{n+1}}{h_2^2}-S_{i,j}^n)\frac{1}{\frac{2}{h_1^2}+\frac{2}{h_2^2}} \end{aligned} \] #### 抛物型方程的求解 ##### 什么是抛物型方程 传导问题的扩散方程是典型的抛物型方程。 方程的形式一般为: \[ \frac{\partial\Phi}{\partial t}=\frac{\partial^2\Phi}{\partial x^2}+S(x,t) \] ##### 显式差分法 下面的式子中的下角标代表着坐标空间(实际上就是坐标),上角标代表着时间空间(代表着时间)。 \[ \begin{aligned} \frac{\Phi_i^{n+1}-\Phi_{i}^{n}}{\Delta t}=\frac{\Phi^n_{i+1}+\Phi^{n}_{i-1}-2\Phi_i^n}{h^2}+S^n \end{aligned} \] \[ \Phi_i^{n+1}=\Phi_i^n+\frac{\Delta t}{h^2}(\Phi^n_{i+1}+\Phi_{i-1}^n-\Phi_i^n)+\Delta tS_i^n \] 有时候可以用一个简写方法: \[ (\delta^2\Phi^n)_i=\Phi^n_{i+1}+\Phi^{n}_{i-1}-2\Phi_i^n \] ##### 隐式差分法 \[ \begin{aligned} \frac{\Phi_i^{n+1}-\Phi_{i}^{n}}{\Delta t}=\frac{\Phi^{n+1}_{i+1}+\Phi^{n+1}_{i-1}-2\Phi_i^{n+1}}{h^2}+S^n_i \end{aligned} \] 化简得到方程:( \(r=\frac{\Delta t}{h^2}\) ) \[ -r\Phi_{i-1}^{n+1}+(1+2r)\Phi_i^{n+1}-r\Phi_{i+1}^{n+1}=\Phi_i^n+S_i^n\Delta t \] 这样,如果知道了n对应的时间的分布(不同i对应的函数值), 那么可以列线性方程组求解下一时刻的函数值。 这个好像只适用于给定了第一类边界条件的情况。就是给定边界的函数值。 如果给定了边界的值，假如说x方向分为了N+1个格点, 第一个和最后一个格点的值已经知道了,那么可以列出一组线性方程。 \[ \begin{bmatrix}1+2r & -r \\-r & 1+2r&-r\\ ...&...&...&...&...&...\\...&...&...&...&...&...\\...&...&...&...&-r&1+2r\\ \end{bmatrix} \begin{bmatrix} \Phi_1^{n+1}\\\Phi_2^{n+1}\\ ..\\.. \\\Phi_{N-1}^{n+1} \end{bmatrix} = \begin{bmatrix} \Phi^n_{1}+S_{1}^n\Delta t+r\Phi_0^{n+1}\\ \Phi^n_{2}+S_{2}^n\Delta t \\ \\ \\ \Phi^n_{N-1}+S_{N-1}^n\Delta t+r\Phi_N^{n+1} \end{bmatrix} \]
平均隐式差分法Crank-Nicolson
类似于隐士差分法,将输运方程的空间二阶导项写为新旧时刻所得值的结合。 \[ \begin{aligned} \frac{\Phi_i^{n+1}-\Phi_{i}^{n}}{\Delta t}=(\frac{1}{2}(\Phi^{n+1}_{i+1}+\Phi^{n+1}_{i-1}-2\Phi_i^{n+1})+\frac{1}{2}(\Phi^{n}_{i+1}+\Phi^{n}_{i-1}-2\Phi_i^{n}))\frac{1}{h^2}+S^n_i \end{aligned} \] 类似于隐式差分,得到一组线性方程 \[ \begin{aligned} &\begin{bmatrix}1+r & -r/2 \\-r/2 & 1+r&-r/2\\ ...&...&...&...&...&...\\...&...&...&...&...&...\\...&...&...&...&-r/2&1+r\\ \end{bmatrix} \begin{bmatrix} \Phi_1^{n+1}\\\Phi_2^{n+1}\\ ..\\.. \\\Phi_{N-1}^{n+1} \end{bmatrix} \\&\quad\quad= \begin{bmatrix} \frac{1}{2}r\Phi_0^n+(1-r)\Phi^n_{1}+\frac{1}{2}r\Phi_2^n+S_{1}^n\Delta t+\frac{r}{2}\Phi_0^{n+1}\\ \frac{1}{2}r\Phi_1^n+(1-r)\Phi^n_{2}+\frac{1}{2}r\Phi_3^n+S_{2}^n\Delta t \\ \\ \\ \frac{1}{2}r\Phi_{N-2}^n+(1-r)\Phi^n_{N-1}+\frac{1}{2}r\Phi_N^n+S_{N-1}^n\Delta t+\frac{r}{2}\Phi_N^{n+1} \end{bmatrix} \end{aligned} \]
算符表示的平均隐式差分法(Crank-Nicolson)方法求解schrodinger方程
我们的schrodinger方程 \[ i\hbar\frac{\partial }{\partial t}\Phi=-\frac{\hbar^2}{2m}\nabla^2\Phi+V\Phi \] 如果取常数\(\hbar=1,2m=1\) 定义算符 \[ H\Phi_i^n=-\frac{1}{h^2}(\Phi_{i+1}^n+\Phi_{i-1}^n-2\Phi_i^n)+V_i\Phi_i^n \] Schrodinger方程可以写为(用了隐式差分法) \[ \frac{\Phi_i^{n+1}-\Phi_i^{n}}{\Delta t}=-i\frac{H}{2}\Phi^{n+1}_{i}-i\frac{H}{2}\Phi^n_{i} \] 上面的式子可以写为 \[ \begin{aligned} \Phi_i^{n+1}&=\frac{1-\frac{i}{2}H\Delta t}{1+\frac{i}{2}H\Delta t}\Phi_i^{n}\\ \Phi_i^{n+1}&=(\frac{2}{1+\frac{i}{2}H\Delta t}-1)\Phi_i^n \end{aligned} \] 如果再定义一个辅助变量: \[ \chi_i^n=\frac{2}{1+\frac{i}{2}H\Delta t}\Phi_i^n \] 则有方程 \[ \begin{aligned} \chi_{i-1}^n+(-2+i\frac{2h^2}{\Delta t}-h^2V_i)\chi_i^n+\chi_{i+1}^n=i\frac{4h^2}{\Delta t}\Phi_i^n \\ \Phi_{i}^{n+1}=\chi_i^n-\Phi_i^n \end{aligned} \] 辅助变量的边界条件定义为\(\chi_0=\chi_N=0\) 。 辅助变量对应的方程组。 \[ \begin{bmatrix}r_1 & 1 \\1 & r_2&1\\ ...&...&...&...&...&...\\...&...&...&...&...&...\\...&...&...&...&1&r_{N-1}\\ \end{bmatrix} \begin{bmatrix} \chi_1^{n}\\\chi_2^{n}\\ ..\\.. \\\chi_{N-1}^{n} \end{bmatrix} = \begin{bmatrix} i\frac{4h^2}{\Delta t}\Phi_1^n-\chi_0^n\\ i\frac{4h^2}{\Delta t}\Phi_2^n\\..\\..\\ i\frac{4h^2}{\Delta t}\Phi_{N-1}^n-\chi_N^n \end{bmatrix} \] ### 第六章分子动力学 #### Verlet算法 算是分子动力学中的最核心的算法了吧。 有两种: * 速度verlet算法 * L.verlet算法 ##### L.verlet算法 牛顿第二定理 \[ \frac{d^2r}{dt^2}=\frac{F(t)}{m} \] 利用三点法改写等式左端。同时速度用前后差分法来确定。 \[ \begin{aligned} \frac{r(t+h)+r(t-h)-2r(t)}{h^2}&=F(t)/m\\ r(t+h)&=2r(t)-r(t-h)+h^2\frac{F(t)}{m}+O(h^4)\\ P(t)&=\frac{m}{2h}(r(t+h)-r(t-h))+O(h^2) \end{aligned} \] 这种算法需要两个启动点,速度方面不需要启动点。 ##### 速度verlet算法 这个是最常用的verlet 算法。 利用匀加速运动的位移公式 \[ \begin{aligned} r(t+h)=r(t)+hv(t)+\frac{h^2}{2}\frac{F(t)}{m}\\ r(t-h)=r(t)-hv(t)+\frac{h^2}{2}\frac{F(t)}{m} \end{aligned} \] 将第二个式子的时间向前推进h。 \[ r(t)=r(t+h)-hv(t+h)+\frac{h^2}{2}\frac{F(t+h)}{m} \] 这个式子和第一个式子相加,可以得到速度的递推公式, 同时坐标的地推公式就用简单的匀加速运动的公式。 \[ \begin{aligned} r(t+h)= r(t)+hv(t)+\frac{h^2}{2}\frac{F(t)}{m}\\ v(t+h)=v(t)+h\frac{F(t)+F(t+h)}{2m} \end{aligned} \] 上面的两个公式就是最常用的速度verlet算法。 #### 边界条件 一般采用的是周期性边界条件,如果A是任意客观测量则周期边界条件的表达式为: \[ A(\vec{x})=A(\vec{x}+\vec{n}L) \] #### 模拟流程 主要就是一点....对温度定义 \[ T= \frac{\bar{E_k}}{\frac{d}{2}Nk_B} \]
蒙特卡洛方法
蒙特卡洛方法算积分, 重要抽样法
考察积分,按照蒙特卡洛方法的思想,我们想随机抽取x，然后得到被积函数的抽样平均值认为这就是被积函数在这个区间上的平均值。 由于f(x)的函数值随着x的变化起伏比较大, 会使得计算误差大。 \[ I = \int_0^xf(x')dx' \] 做变量代换 后的积分可以写为 \[ I = \int_0^x\frac{f(x')}{w(x')}w(x')dx' \] 如果取\(y(x)=\int_0^xw(x')dx'\),并且要求\(\int_0^xw(x')dx'=1\) \[ I = \int_0^1\frac{f(x'(y))}{w(x'(y))}dy \] 至此相当于对积分做了一次积分变换。 然后我们在\([0,1]\)的区间抽取N个样本y, 把样本叫做\(y_i,y_i\in[0,1]\)。则积分可以表示为 \[ I = \frac{1}{N} \Sigma_1^N \frac{f(x'(y_i))}{w(x'(y_i))} \]
这个积分方法就叫做重要抽样法。 解释这个方法为什么是对的 ,这段话是定性的。 当N很大是,可以看做\(y_i\)在\([0,1]\)上面均匀分布。在区间\(y\to y+dy\)内, \(y_i\)落入其中的概率(的期望)为\(\frac{dy}{1}\), 也就是说, 抽到应变量为\(\frac{f(x'(y))}{w(x'(y))}\),(严格来讲应该是\(\frac{f(x'(y+\epsilon dy))}{w(x'(y+\epsilon dy))}\))的概率(的期望)为\(dy\)。(求和化为积分的具体过程)那么,遍历整个积分区间\([0,1]\)。积分是抽取到应变量的平均值也是积分值。
重要抽样法的数值方法1,差分方法
当我们选取的函数\(w(x)\)不方便积分时, 重要抽样法就难以进行。 在\([0,1]\)的区间中选取\(M+1\)个y点。可以表示为\(y^{i}=\frac{i}{M},i\in\{0,1...,M\}\)。 由于 \[ \frac{dy}{dx} = w(x) \] 化简为差分方程 \[ \frac{y^{i+1}-y^{i}}{x^{i+1}-x^i}=w(x^i) \] \[ x^{i+1}=x^i+\frac{1}{Mw(x^i)} \] 利用上面的递推关系,结合起始点\(x^0=0\)，可以得到\(x^i\)的具体值。 当然,选取w(x)时, 这里默认是保证了\(\int_0^xw(x')dx'=1\)
蒙特卡洛方法的数值方法2, Von_Neumann 舍选法
我们想让抽取的x样本随x的分布按照\(w(x)\)进行。选取函数\(w'(x)\),满足\(w'(x)>w(x)\)。 这个抽样这样进行:
- 在区间\([0,x]\)中抽取随机数\(x_i\)
- 计算\(w'(x_i)\)和\(w(x_i)\)。
- 在区间\([0,1]\)之间生成随机数\(\eta\)
- 比较\(\frac{w(x_i)}{w'(x_i)}\)和\(\eta\)
- 若\(\frac{w(x_i)}{w'(x_i)}>\eta\)则这个\(x_i\)保留, 否则再抽取下一个\(x_i\) 实际操作时候\(w'(x)=max_{x'\in [0,x]} w(x')\) 是一个常数。 按照这个流程,保留下来的\(x_i\)都倾向于集中在\(w(x_i)\)大的地方。 当然,选取w(x)时, 这里默认是保证了\(\int_0^xw(x')dx'=1\)
重要抽样法计算多维积分, metropolis方法
在多维空间中,这样选取随机点。 1. 确定空间中的重要性函数w(X)。 2.
确定起始点\(X_0\),
这个是一个列向量,或者说是一个坐标 3. 下一步\(X_{i+1}=X_i+\delta_r\)。 其中\(\delta_r\)是指对于X的每一个坐标都变化一个在\([-\delta,\delta]\)中的随机数。 4. 在\([0,1]\)中选取一个随机数 \(\eta\), 如果\(\frac{w(X_{i+1})}{w(X_i)}>\eta\)则保留\(X_{i+1}\)。否则\(X_{i+1}=X_i\)
按照这个方法选取的随机数会倾向于到\(w(X)\)大的地方去。 当然,选取w(x)时,
这里默认是保证了\(\int_0^xw(x')dx'=1\) *
一个metropolis方法的例子(matlab) 是老师写的程序,
这个程序考察了规定重要性\(w(x)=\frac{6}{5}(1-\frac{1}{2}x^2)\)
。第一个点取的是\(x=0.9\),
然后将随机取得点的分布画出来和\(w(x)\)进行对比。
但是实际上这个程序没有体现出metropolis方法中的"行走"。
fxt的代码是这样处理行走中的边界问题的。(见下一段代码),对游走的点分了类，边界的点有相应的游走步长。
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
53function metropolis_2
clear all
close all
clc
x1=0.9;
delta=5.5;
w=@(x1) 6/5*(1-0.5*x1.^2); %概率密度函数，与例题8相同
n=0;
N=5e5;
for i=1:N
[x1]=my_metropolis(x1,delta,w);
n=n+1;
xx(n)=x1; %将经历的x1值存储起来
end
%下面的程序判断metropolis算法生成的随机数分布是否符合分布函数
%========================================================
Num=201; %区间大小
Max=max(abs(xx))+10^-9 %边界
Min=min(xx)*1.0 %边界
y=zeros(1,Num);
x=linspace(Min,Max,Num);
h=x(2)-x(1);
for i=1:length(xx)
index=fix((xx(i)-Min)/h)+1; %判定xx(i)处于哪个子区间，fix为舍去小数取整运算
y(index)=y(index)+1; %该子区间的随机数数目+1
end
f=@(x1) 6/5*(1-0.5*x1.^2); %精确值绘图
figure,fplot(f,[Min,Max],'r');
hold on;
S=interg(Min,Max,f)
y=y/length(xx)/h*S; %生成的随机数的数目为n(length(xx))，走的总步数为N
plot(x(1:Num-1),y(1:Num-1),'-*');title('验证随机数的分布');
legend('精确分布','随机数的分布')
hold off
function S=interg(a,b,w) %w函数在[a,b]上的积分
n=10000;
h=(b-a)/n;
xx_h=a:h:b;
ww=w(xx_h);
%S=sum(w(xx_h).*(b-a)/n)-w(a)*(b-a)/2/n-w(b)*(b-a)/2/n; %梯形算法
S=(4*sum(ww(2:2:n))+2*sum(ww(3:2:n-1))+ww(1)+ww(n+1))*h/3; %simpson积分公式
%=========================================================================
function [x1] = my_metropolis(x1,delta,w)
x1t=rand; %[0,1]区间产生随机数
r=w(x1t)/w(x1);
if rand<r;
x1=x1t;
end
1 | function [xnew,typenew] = MH(x,delta,type,w) |
蒙特卡洛方法计算Ising 模型
- 简单介绍Ising模型 二维Ising模型就是一个\(N_x\times N_y=N\)的格点每一个格点上面都有一个自旋, 记为\(S_{i,j}\) 系统有一个哈密顿量可以写为\(H = -J\Sigma_{<\alpha \beta>}S_{\alpha}S_{\beta}-B\Sigma_{\alpha}S_{\alpha}\) 上面的式子中,\(<\alpha \beta>\)指的是对所有相邻的连线求和。当然就此可以引入两种边界条件 (周期边界条件, 螺旋边界条件) 一个分布的权重(重要性)\(w(S)=\frac{e^{-H(S)}}{Z}\),其中\(Z = \Sigma e^{-H(S)}\) 。 二维Ising 模型经过计算可以得到物理量和统计学期望之间的关系。 \[ \begin{aligned} M &= <\Sigma_{\alpha = 1}^N S_{\alpha}>\\ \chi &= <\Sigma_{\alpha = 1}^N S_{\alpha}^2>-<\Sigma_{\alpha = 1}^N S_{\alpha}>^2\\ E &= <H(S)>\\ C_B &= <H^2>-<H>^2 \end{aligned} \]
Ising 模型用metropolis方法的具体操作 1. 选择一个初始的S的分布(构型) 2. 选择一个格点{i,j}。 3. 改变{i,j}上的自旋(使他反号) 4. 计算改变后的系统能量\(H'=H+(2Jf+B)S_{\alpha}\) ,则改变前后的的权重比值可以写为\(r=\frac{w(S')}{w(S)}=e^{-S_{\alpha}(2Jf+B)}\), 其中\(f = (S_{i-1,j}+S_{i+1,j}+S_{i,j+1},S_{i,j-1})\) 。实际上f只能有5个可以取得值(从-4到4,间隔2取值),然后\(S_{\alpha}\)也只有两个可以取的值，综合来看, r只有十个值，将这十个值提前存储在列表里面可以节省计算时间。 5. 从\([0,1]\)区间里面取值\(\eta\),如果\(\eta<r\),则认为变化后的分布\(S_{i+1}\)就是\(S_i\)翻转了一个自旋；如果\(\eta> r\), 则认为\(S_{i+1}\)就是\(S_i\)。 6. 在进行了很多步之后就可以开始计算物理量了。
- 下面是Ising模型的例子,
选择格子在水平和竖直方向上的大小都为50，耦合强度和外磁场强度分别为J=0.3，
B=0.热化扫描50次（此时不计算可观测量），数据分成10个小组，每组50个样本，每个样本的抽样间隔为5.
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
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170function Exapmple_Ising_model
clear all
close all
clc
NX=6;NY=6; %格子尺寸
J=[0.0:0.05:0.3,0.305:0.005:0.605,0.61:0.05:0.92];
%耦合强度，并不等间隔的原因是为了加强0.3-0.6段的计算
B=0.0; %外磁场强度
N_therm=50; %热化扫描步数
N_group=10; %数据分组数目
N_size=50; %每组样本数
N_freq=5; %抽样间隔步数
S=zeros(NX,NY);
[N,S]=initialize(NX,NY); %初始化格子的自旋
S1=S
for iii=1:length(J)
JJ=J(iii); %依次取不同的J值,即不同的耦合强度
R=zeros(5,2); %R数组用于存放接受概率值，有10个矩阵元，即PPT中的(2)式
R=flip_probs(JJ,B);% %初始化R数组
for i=1:N_therm %做热化扫描
[S,accept]=thermal_sweep(NX,NY,S,R);
%accept/N %反转操作接受率
end
S2=S;
Sum_M=0; Sum_M2=0; Sum_SigM=0; %初始化待求物理量M, E，Chi，CB
Sum_E=0; Sum_E2=0; Sum_SigE=0;
Sum_Chi=0; Sum_Chi2=0;
Sum_CB=0; Sum_CB2=0;
for igroup=1:N_group
Group_M=0; Group_M2=0;
Group_E=0; Group_E2=0;
for sweep=1:N_size*N_freq %完成此轮循环，就完成了一组数值的计算
[S,accept]=thermal_sweep(NX,NY,S,R); %抽样前，先做热化扫描
%accept/N %反转操作接受率
if mod(sweep,N_freq)==0
[Magic,Energy]=Magic_Energy(B,JJ,NX,NY,S); %求出一种情况下的磁化强度和能量E
%一次抽样，获得M和E
Group_M=Group_M+Magic;
Group_M2=Group_M2+Magic^2;
Group_E=Group_E+Energy;
Group_E2=Group_E2+Energy^2;
end
end
%===================================
Group_M=Group_M/N_size; %组内抽样的平均，即期望值，PPT中的（1）式，教材中的(8.21a)
Group_M2=Group_M2/N_size;
Group_E=Group_E/N_size;
Group_E2=Group_E2/N_size;
Chi=Group_M2-Group_M^2; %(8.21b),一个组里面求得的磁化率
CB=Group_E2-Group_E^2; %(8.21d)
Sum_M=Sum_M+Group_M;
Sum_E=Sum_E+Group_E;
Sum_Chi=Sum_Chi+Chi;
Sum_CB=Sum_CB+CB;
end
E=Sum_E/N_group;
M=Sum_M/N_group;
Chi=Sum_Chi/N_group; %几个组求得的磁化率的平均
CB=Sum_CB/N_group;
CBB(iii)=CB; %用于存储不同J值对应的CB比热容
EE(iii)=E/N; %N为总格子数
MM(iii)=M/N;
Chi_1(iii)=Chi/N;
end %对应于for iii=1:length(J)
figure;plot(1./J,CBB,'r*-') %注意：由于1/0为无穷大，程序没有画出该点
title('比热随耦合强度的倒数的变化关系')
xlabel('1/J'); ylabel('C_B');
figure;plot(1./J,EE,'b') %J值大，意味着温度低，M值大
title('能量随耦合强度的倒数的变化关系')
xlabel('1/J'); ylabel('能量E');
figure;plot(1./J,MM,'r') %注意，磁化强度可以朝不同方向，因此可以为负
title('磁化强度随耦合强度的倒数的变化关系')
xlabel('1/J'); ylabel('磁化强度M');
figure;plot(1./J,Chi_1,'k')
title('磁化率随耦合强度的倒数的变化关系')
xlabel('1/J'); ylabel('磁化率Chi');
S
function [N,S]=initialize(NX,NY) %初始化每个格子的自旋
N=NX*NY;
for i=1:NY
for j=1:NX
if rand < 0.5
S(i,j)=1;
else
S(i,j)=-1;
end
end
end
function [Magic,Energy]=Magic_Energy(B,JJ,NX,NY,S) %求出一种情况下的磁化强度和能量E
Magic=0;
Sum_ss=0;
for i=1:NY
if i>1
IM=i-1;
else
IM=NY;
end
%====================
for j=1:NX
if j>1
JM=j-1;
else
JM=NX;
end
%==================
Magic=Magic+S(i,j);
Sum_ss=Sum_ss+S(i,j)*(S(IM,j)+S(i,JM));
%注意上式求和实际上是进行(8.18)式第一项的运算，但S(i,j)只与左格点和上格点相乘求和
%是因为如果上下左右四个格点都相乘求和，则存在重复计算的现象
end
end
Energy=-JJ*Sum_ss-B*Magic; %(8.18)式
function R=flip_probs(JJ,B) %设定翻转判定矩阵
for i=1:5
R(i,2)=exp(-2*(JJ*(2*i-6)+B)); %Sa为1
%自旋反转判定因子，page 147页的（8.23）式，PPT中的（2）式，R(:,2)代表中心点自旋向上
R(i,1)=exp(2*(JJ*(2*i-6)+B)); %分别对应于f=-4，-2, 0, 2, 4
end
function [S,accept]=thermal_sweep(NX,NY,S,R)
%热化扫描函数，每次热化扫描后，所有的点都被遍历一次
accept=0;
for i=1:NY
if i<NY %采用周期性边界条件，此if循环给定(i,j)粒子下面的相邻格点（i+1,j）
IP=i+1;
else
IP=1;
end
if i>1 %采用周期性边界条件，此if循环给定(i,j)粒子上面的相邻格点（i-1,j）
IM=i-1;
else
IM=NY;
end
%==================
for j=1:NX
if j<NX %采用周期性边界条件，此if循环给定(i,j)粒子右边的相邻格点（i,j+1）
JP=j+1;
else
JP=1;
end
if j>1 %采用周期性边界条件，此if循环给定(i,j)粒子左边的相邻格点（i,j-1）
JM=j-1;
else
JM=NX;
end
%==================
spin=S(i,j);
f=S(IP,j)+S(IM,j)+S(i,JP)+S(i,JM);
if rand>R(3+f/2,(3+spin)/2) %如果反转概率比随机数小，则不反转
%注意：f的取值只有-4，-2,0,2,4五种，spin只有-1,1两种，
%因此用R（5,2）二维数组即可以判定是否反转
else
S(i,j)=-spin; %反转(i,j)点的自旋
accept=accept+1;
end
end
end
典型的例子
在作业文件夹中有一个"自定义函数"的子文件夹,里面有很多实用的, 暂时没有搬到笔记里面来。
4-RK算法
写了一个函数文件,可以以后调用。 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%这里是想写一个给定初态,分段,分段数目,函数迭代表达式,就可以返回函数自变量和因变量。
%f,可以是多因变量的,但是目前只支持单自变量。f的定义格式可以是f = @(x,y)[dy1/dx(x,y) ,dy2/dx(x,y) ...]
%y 是因变量,这里要写成行变量的形式。[y1,y2...]
%x_start是x的起始点
%x_end是x的终止点
%n 是分段的数目,相当于自变量列表的数目是n+1
%支持倒序解
function [x,y]= RK_solution(x_start,x_end,y_start,n,f)
h=(x_end-x_start)/n;
x=x_start:h:x_end;
%通过这个变量(length)得到函数的因变量。
length = size(y_start,2);
%将因变量定义为元素为0的矩阵。
y = zeros(n+1,length);
%将函数因变量的第一行写为初态。
y(1,:)=y_start;
%利用R-K算法迭代
for i = 1:1:n
xn=x(i);
yn=y(i,:);
K1 = f(xn,yn);
K2 = f(xn+h/2,yn+h/2.*K1);
K3= f(xn+h/2,yn+h/2.*K2);
K4= f(xn+h,yn+h.*K3);
y(i+1,:)=yn+h/6.*(K1+2.*K2+2.*K3+K4);
end
end
Simpson38积分法
1 | % 定义simpson38算法的积分,算法可以做到给定函数和自变量数组就可以,需要注意的是x中的元素个数一定要是。3的倍数+1 |
简单搜索发结合二分法
1 | % 定义一个二分法结合简单搜索发寻找哦函数f从xstart到xend之间所有的根,并且寻找到的根的数目可以自定义为nmax. |
薛定谔方程的本征值问题
[[计算物理#求解薛定谔方程的本征值问题|前面章节讲到的打靶法]]
1 | %[energy,x,y] = find_energy(x_min,x_max,h,v,gamma,n,delta,delta_judge,epsilon,e_start,e_end,e_max,y1a,y1b) |
索末菲方法
1 | %sommerfeld_energy = sommerfeld(n,v,gamma,h,delta,epsilon,x_l,x_r,e_start) |