四阶龙格库塔法解一维扩散方程

龙格库塔法是解ode的利器，最近有个project是用龙格库塔法解pde，做的时候网上这方面能搜到的资料还不多，传一下自己的解法。

核心思路

直线法（method of line）把pde化为ode，离散化为空间n个格点上的一元函数，用向量表示
u(x,ti)→{u1(t),u2(t)...,un(t)}iT=u⃗iu(x,t_i)→\{u_1(t),u_2(t)...,u_n(t)\}_i^T=\vec{u}_i u(x,ti)→{u1(t),u2(t)...,un(t)}iT=ui
空间的二阶导用中心差分格式处理，写成矩阵形式
∂u2∂2x=1Δx2(−2100⋯01−210⋯⋮01−21⋯⋮⋮⋮⋮⋮⋯1000⋯1−2)⋅{u1(t)u2(t)⋮⋮un(t)}=∂u∂t=C⋅u⃗i\frac{\partial u^2}{\partial^2 x}=\frac{1}{\Delta x^2} \begin{pmatrix} -2 & 1 & 0 & 0 &\cdots & 0 \\ 1 & -2 & 1 & 0 &\cdots & \vdots \\ 0 & 1 & -2 & 1 & \cdots & \vdots \\ \vdots & \vdots & \vdots & \vdots & \cdots & 1 \\ 0 & 0 & 0 & \cdots & 1 & -2 \end{pmatrix}·\left \{ \begin{matrix} u_1(t)\\u_2(t)\\\vdots \\\vdots \\u_n(t) \end{matrix} \right\} =\frac{\partial u}{\partial t} =C·\vec{u}_i ∂2x∂u2=Δx21⎝⎜⎜⎜⎜⎜⎜⎜⎛−210⋮01−21⋮001−2⋮0001⋮⋯⋯⋯⋯⋯10⋮⋮1−2⎠⎟⎟⎟⎟⎟⎟⎟⎞⋅⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧u1(t)u2(t)⋮⋮un(t)⎭⎪⎪⎪⎪⎪⎪⎬⎪⎪⎪⎪⎪⎪⎫=∂t∂u=C⋅ui
扩散方程的时间一阶偏导用4阶龙格库塔方法转化为线性方程
四阶龙格库塔常用格式

∂u∂t=∂2u∂x2→ui+1−uih=16(k1+k2+k3+k4)=C⋅ui⃗\frac{\partial u}{\partial t} = \frac{\partial^2u}{\partial x^2}\to \frac{u_{i+1}-u_{i}}{h}=\frac{1}{6}(k_1+k_2+k_3+k_4) =C·\vec{u_i} ∂t∂u=∂x2∂2u→hui+1−ui=61(k1+k2+k3+k4)=C⋅ui
4.中间量取值
这是最关键的一步。龙格库塔解ode之所以方便，因为f是表达式已知的显函数，直接代入完事。但pde里面偏导表达式如果有还好办，没有的话像扩散方程得换条思路。注意到k是时间偏导的函数值，取决于x和y。在扩散方程中k其实是空间二阶偏导的在某一点的函数值。虽然定义上某一点y的二阶偏导由x,y共同决定，但实际上把二阶偏导视为微分算符的话，偏导的值由y便可确定。换言之，扩散方程中u是一个隐函数，不同的k对应不同的u。\begin{matrix} {\color{Red} 这是最关键的一步。}\\ {\color{Red} 龙格库塔解ode之所以方便，因为f是表达式已知的显函数，直接代入完事。}\\ {\color{Red}但pde里面偏导表达式如果有还好办，没有的话像扩散方程得换条思路。} \\ {\color{Red}注意到k是时间偏导的函数值，取决于x和y。} \\ {\color{Red}在扩散方程中k其实是空间二阶偏导的在某一点的函数值。} \\ {\color{Red}虽然定义上某一点y的二阶偏导由x,y共同决定，}\\ {\color{Red}但实际上把二阶偏导视为微分算符的话，偏导的值由y便可确定。}\\ {\color{Red}换言之，扩散方程中u是一个隐函数，不同的k对应不同的u。}\\ \end{matrix} 这是最关键的一步。龙格库塔解ode之所以方便，因为f是表达式已知的显函数，直接代入完事。但pde里面偏导表达式如果有还好办，没有的话像扩散方程得换条思路。注意到k是时间偏导的函数值，取决于x和y。在扩散方程中k其实是空间二阶偏导的在某一点的函数值。虽然定义上某一点y的二阶偏导由x,y共同决定，但实际上把二阶偏导视为微分算符的话，偏导的值由y便可确定。换言之，扩散方程中u是一个隐函数，不同的k对应不同的u。

k1=C⋅ui⃗,k2=C⋅（ui⃗+hk12）,k3=C⋅（ui⃗+hk22）,k4=C⋅（ui⃗+hk32）k_1=C·\vec{u_i},k_2=C·（\vec{u_i}+\frac{hk_1}{2} ）,k_3=C·（\vec{u_i}+\frac{hk_2}{2} ）,k_4=C·（\vec{u_i}+\frac{hk_3}{2} ） k1=C⋅ui,k2=C⋅（ui+2hk1）,k3=C⋅（ui+2hk2）,k4=C⋅（ui+2hk3）

迭代
至此，我们得到一个在时间方向上迭代的前向差分格式
ui+1=ui+h6(k1+k2+k3+k4){u_{i+1}=u_{i}}+\frac{h}{6}(k_1+k_2+k_3+k_4) ui+1=ui+6h(k1+k2+k3+k4)
C++用一个循环实现，调用eigen线代库进行矩阵运算，vector容器定义循环的列向量，ofstream把每次循环结果列向量的值装进excel表。

计算结果

边界与初始条件，以及计算范围与步长
D=0.1u(0,t)=1013,u(x>0,0)=0t∈[0,5],tsteps=10000;x∈[0,1],xsteps=50D=0.1\\ u(0,t)=10^{13},u(x>0,0)=0 \\ t∈[0,5],tsteps=10000;x∈[0,1],xsteps=50 D=0.1u(0,t)=1013,u(x>0,0)=0t∈[0,5],tsteps=10000;x∈[0,1],xsteps=50
小步长多步计算，基本等于恒定源扩散的解析解——余误差函数

缩小范围，减少步数，计算误差
D=0.01t∈[0,3],tsteps=300;x∈[0,1],xsteps=30D=0.01\\ t∈[0,3],tsteps=300;x∈[0,1],xsteps=30 D=0.01t∈[0,3],tsteps=300;x∈[0,1],xsteps=30
龙格库塔算出来的误差随函数值增大而增大，空间越往后误差越大，时间越往后误差越小

对比与欧拉法的误差

误差龙格库塔-欧拉<0，说明确实比欧拉法精确。