线性抛物方程的通用有限差分格式及代码（未完）

题目

通用差分格式

现在是type1 dirichlet边界条件
将显示格式和隐式格式写成一个格式，当$\theta =1$是隐式格式，$\theta =0$是显式格式，$\theta =1/2$是CN格式。

离散区域 $M=(L-l)/\Delta x,N=T/\Delta t,k=1...M+1$,n=1…N+1
离散初值 $u_k^1=g(k\ast \Delta x),k=1,...,M+1$
离散时间项 $$\frac{u_k^{(n+1)}-u_k^{(n)}}{\Delta t}=u_t(k,n)，n=1,...,N$$
离散空间项
隐式格式$$\frac{{\delta }^2{u}_k^{n+1}}{\Delta x^2}=u_{xx}(k,n)，K=2,...,M$$
$$\frac{\delta u_k^{n+1}}{2\Delta x}=u_x(k,n)，K=2,...,M$$
显式格式$$\frac{{\delta }^2{u}_k^n}{\Delta x^2}=u_{xx}(k,n)，n=1,...,N$$ $$\frac{\delta u_k^n}{2\Delta x}=u_x(k,n)，K=2,...,M$$
将差分形式代入
$$\frac{u_k^{(n+1)}-u_k^{(n)}}{\Delta t}-\theta \nu \frac{{\delta }^2{u}_k^{n+1}}{\Delta x^2}-(1-\theta )\frac{{\delta }^2{u}_k^n}{\Delta x^2}+\theta a\frac{\delta u_k^{n+1}}{2\Delta x}+(1-\theta )a\frac{\delta u_k^{}}{2\Delta x}+b\theta u_k^{n+1}+b(1-\theta )u_k^n=\theta F_k^{n+1}+(1-\theta )F_k^n$$

$$u_k^{(n+1)}-u_k^{(n)}-\Delta t\theta \nu \frac{{\delta }^2{u}_k^{n+1}}{\Delta x^2}-\Delta t(1-\theta )\frac{{\delta }^2{u}_k^n}{\Delta x^2}+\Delta t\theta a\frac{\delta u_k^{n+1}}{2\Delta x}+\Delta t(1-\theta )a\frac{\delta u_k^{}}{2\Delta x}+\Delta tb\theta u_k^{n+1}+\Delta tb(1-\theta )u_k^n=\Delta t\theta F_k^{n+1}+\Delta t(1-\theta )F_k^n$$
令$r=\nu \frac{\Delta t}{\Delta x^2},R=a\frac{\Delta t}{2\Delta x}$移项
$$u_k^{n+1}-r\theta {\delta }^2{u}_k^{n+1}+R\theta \delta u_k^{n+1}+\Delta tb\theta u_k^{n+1}=u_k^n+r(1-\theta ){\delta }^2{u}_k^n-R(1-\theta )\delta u_k^n-\Delta tb(1-\theta )u_k^n+\Delta t\theta F_k^{n+1}+\Delta t(1-\theta )F_k^n$$
写成矩阵 $Q_1{u}_k^{n+1}=Q_2{u}_k^n+F$
$$A_1=\left(\begin{array}{cccc}-2& 1& 0& \cdots \\ 1& -2& 1& \cdots \\ ⋮& ⋮& \ddots & \end{array}\right)$$
$$A_2=\left(\begin{array}{cccc}0& 1& 0& \cdots \\ -1& 0& 0& \cdots \\ ⋮& ⋮& \ddots & \end{array}\right)$$
$$Q_1=(1+\Delta tb\theta )E-r\theta A_1+R\theta A_2$$
$$Q_2=(1-\Delta tb(1-\theta ))E+r(1-\theta )A_1-R(1-\theta )A_2$$
$$F=\left(\begin{array}{c}r\theta u_1^{n+1}+R\theta u_1^{n+1}+r(1-\theta )u_1^n+R(1-\theta )u_1^n\\ 0\\ ⋮\\ 0\\ r\theta u_{M+1}^{n+1}+R\theta u_{M+1}^{n+1}+r(1-\theta )u_{M+1}^n+R(1-\theta )u_{M+1}^n\end{array}\right)+\Delta t\theta F_k^{n+1}+\Delta t(1-\theta )F_k^n$$

代码