当前位置: 首页 > news >正文

偏微分方程算法之混合边界条件下的差分法

目录

一、研究目标

二、理论推导

三、算例实现

四、结论


一、研究目标

        我们在前几节中介绍了Poisson方程的边值问题,接下来对椭圆型偏微分方程的混合边值问题进行探讨,研究对象为:

\left\{\begin{matrix} -(\frac{\partial^{2}u(x,y)}{\partial x^{2}}+\frac{\partial^{2}u(x,y)}{\partial y^{2}})=f(x,y),(x,y)\in\Omega,\space\space(1)\\ \frac{\partial u(x,y)}{\partial \mathbf{n}}+\lambda(x,y)u=\mu(x,y),(x,y)\in\partial \Omega=\Gamma \end{matrix}\right.

其中,\Omega为矩形区域\Omega=[(x,y)|a\leqslant x\leqslant b,c\leqslant y\leqslant d]f(x,y)\Omega上的连续函数,\mathit{\mathbf{n}}\Gamma上的单位法向量,从而\frac{\partial u(x,y)}{\partial\mathbf{n}}表示方向导数,\lambda(x,y),\mu(x,y)\Gamma的连续函数且\lambda(x,y)非负。对于矩形区域\Omega而言,其边界上的法向量没有统一的表达式,需要对四条边界线段分别讨论。可见:

        在左边界\mathbf{n}_{1}=(-1,0),边界条件为(-\frac{\partial u}{\partial x}+\lambda(x,y)u)|_{(a,y)}=\mu_{1}(a,y)

        在右边界\mathbf{n}_{2}=(1,0),边界条件为(\frac{\partial u}{\partial x}+\lambda(x,y)u)|_{(b,y)}=\mu_{2}(b,y)

        在下边界\mathbf{n}_{3}=(0,-1),边界条件为(-\frac{\partial u}{\partial y}+\lambda(x,y)u)|_{(x,c)}=\mu_{3}(x,c)

        在上边界\mathbf{n}_{4}=(0,1),边界条件为(\frac{\partial u}{\partial y}+\lambda(x,y)u)|_{(x,d)}=\mu_{4}(x,d)

        于是,公式(1)可以具体写为:

\left\{\begin{matrix} -(\frac{\partial^{2}u(x,y)}{\partial x^{2}}+\frac{\partial^{2}u(x,y)}{\partial y^{2}})=f(x,y),(x,y)\in(a,b)\times (c,d),\\ (\frac{\partial u(x,y)}{\partial x}-\lambda(x,y)u)|_{(a,y)}=\varphi_{1}(y),c\leqslant y\leqslant d,\\ (\frac{\partial u(x,y)}{\partial x}+\lambda(x,y)u)|_{(b,y)}=\varphi_{2}(y),c\leqslant y\leqslant d,\\ (\frac{\partial u(x,y)}{\partial y}-\lambda(x,y)u)|_{(x,c)}=\Psi_{1}(x),a\leqslant x\leqslant b,\\ (\frac{\partial u(x,y)}{\partial y})+\lambda(x,y)u)|_{(x,d)},\Psi_{2}(x),a\leqslant x\leqslant b \end{matrix}\right.

二、理论推导

        首先进行矩形区域\Omega的等距剖分,得到各网格节点(x_{i},y_{j}),且x_{i}=a+i\cdot\Delta x,i=0,1,\cdot\cdot\cdot,m\Delta x=(b-a)/my_{j}=c+j\cdot\Delta y,j=0,1,\cdot\cdot\cdot,n\Delta y=(c-d)/n。然后弱化原方程,使其仅在离散节点上成立,从而有

\left\{\begin{matrix} -(\frac{\partial^{2}u(x,y)}{\partial x^{2}}+\frac{\partial^{2}u(x,y)}{\partial y^{2}})|_{(x_{i},y_{j})}=f(x_{i},y_{j}),1\leqslant i\leqslant m-1,1\leqslant j\leqslant n-1,\\ \frac{\partial u(x,y)}{\partial x}|_{(x_{0},y_{j})}-\lambda(x_{0},y_{j})u(x_{0},y_{j})=\varphi_{1}(y_{j}),0\leqslant j\leqslant n, \\ \frac{\partial u(x,y)}{\partial x}|_{(x_{m},y_{j})}+\lambda(x_{m},y_{j})u(x_{m},y_{j})=\varphi_{2}(y_{j}),0\leqslant j\leqslant n, \\ \frac{\partial u(x,y)}{\partial y}|_{(x_{i},y_{0})}-\lambda(x_{i},y_{0})u(x_{i},y_{0})=\Psi_{1}(x_{i}),0\leqslant i\leqslant m, \\ \frac{\partial u(x,y)}{\partial y}|_{(x_{i},y_{n})}+\lambda(x_{i},y_{n})u(x_{i},y_{n})=\Psi_{2}(x_{i}),0\leqslant i\leqslant m. \end{matrix}\right.

        将上式中的一阶、二阶偏导数分别用关于一阶导数的中心差商和关于二阶导数的中心差商来近似,得

\left\{\begin{matrix} -(\frac{u(x_{i-1},y_{j})-2u(x_{i},y_{j})+u(x_{i+1},y_{j})}{\Delta x^{2}}+\frac{u(x_{i},y_{j-1})-2u(x_{i},y_{j})+u(x_{i},y_{j+1})}{\Delta y^{2}})\\ =f(x_{i},y_{j})+O(\Delta x^{2}+\Delta y^{2}),1\leqslant i\leqslant m-1,1\leqslant j\leqslant n-1,\\ \frac{u(x_{1},y_{j})-u(x_{-1},y_{j})}{2\Delta x}-\lambda(x_{0},y_{j})u(x_{0},y_{j})=\varphi_{1}(y_{j})+O(\Delta x^{2}),0\leqslant j\leqslant n,\\ \frac{u(x_{m+1},y_{j})-u(x_{m-1},y_{j})}{2\Delta x}+\lambda(x_{m},y_{j})u(x_{m},y_{j})=\varphi(y_{j})+O(\Delta x^{2}),0\leqslant j\leqslant n,\\ \frac{u(x_{i},y_{1})-u(x_{i},y_{-1})}{2\Delta y}-\lambda(x_{i},y_{0})u(x_{i},y_{0})=\Psi_{1}(x_{i})+O(\Delta y^{2}),0\leqslant i\leqslant m,\\ \frac{u(x_{i},y_{n+1})-u(x_{i},y_{n-1})}{2\Delta y}+\lambda(x_{i},y_{n})u(x_{i},y_{n})=\Psi_{2}(x_{i})+O(\Delta y^{2}),0\leqslant i\leqslant m. \end{matrix}\right.

        用数值解代替精确解并忽略高阶项,得到以下数值格式:

\left\{\begin{matrix} -(\frac{u_{i-1,j}-2u_{i,j}+u_{i+1,j}}{\Delta x^{2}}+\frac{u_{i,j-1}-2u_{i,j}+u_{i,j+1}}{\Delta y^{2}})=f_{i,j},1\leqslant i\leqslant m-1,1\leqslant j\geqslant n-1,\\ u_{1,j}-u_{-1,j}=2\Delta x(\varphi_{1}(y_{j})+\lambda_{0,j}u_{0,j}),0\leqslant j\leqslant n,\\ u_{m+1,j}-u_{m-1,j}=2\Delta x(\varphi_{2}(y_{j})-\lambda_{m,j}u_{m,j}),0\leqslant j\leqslant n,\\ u_{i,1}-u_{i,-1}=2\Delta y(\Psi_{1}(x_{i})+\lambda_{i,0}u_{i,0}),0\leqslant i\leqslant m,\\ u_{i,n+1}-u_{i,n-1}=2\Delta y(\Psi_{2}(x_{i})-\lambda_{i,n}u_{i,n}),0\leqslant i\leqslant m. \end{matrix}\right.(2)

其中,\lambda_{i,j}=\lambda(x_{i},y_{j}),f_{i,j}=f(x_{i},y_{j})。该格式的局部截断误差为O(\Delta x^{2}+\Delta y^{2}),同时需要处理其中的下标越界问题。由于函数的连续性,我们认为公式(2)中的第1式对i=0也成立,即有

-(\frac{u_{-1,j}-2u_{0,j}+u_{1,j}}{\Delta x^{2}}+\frac{u_{0,j-1}-2u_{0,j}+u_{0,j+1}}{\Delta y^{2}})=f_{i,j},1\leqslant j\leqslant n-1

        再成公式(2)的第2式中解出u_{-1,j}=u_{1,j}-2\Delta x(\varphi_{1}(y_{j})+\lambda_{0,j}u_{0,j})代入上式,得

-\frac{2u_{1,j}-2(1+\Delta x\lambda_{0,j})u_{0,j}}{\Delta x^{2}}-\frac{u_{0,j-1}-2u_{0,j}+u_{0,j+1}}{\Delta y^{2}}=f_{i,j}-\frac{2}{\Delta x}\varphi_{1}(y_{j}),1\leqslant j\leqslant n-1

同样,设公式(2)中的第1式分别对i=m,j=0,j=n成立,再分别从公式(2)的第3、4、5式中解出u_{m+1,j},u_{i,-1},u_{i,n+1}代入前面刚得到的方程,就可以处理掉越界下标,得到以下格式:

\left\{\begin{matrix} -\frac{u_{i-1,j}-2u_{i,j}+u_{i+1,j}}{\Delta x^{2}}-\frac{u_{i,j-1}-2u_{i,j}+u_{i,j+1}}{\Delta y^{2}}=f_{i,j},1\leqslant i\leqslant m-1,1\leqslant j\leqslant n-1,\\ -\frac{2u_{1,j}-2(1+\Delta x\lambda_{0,j})u_{0,j}}{\Delta x^{2}}-\frac{u_{0,j-1}-2u_{0,j}+u_{0,j+1}}{\Delta y^{2}}=f_{0,j}-\frac{2}{\Delta x}\varphi_{1}(y_{j}),1\leqslant j\leqslant n-1,\\ -\frac{2u_{m-1,j}-2(1+\Delta x\lambda_{m,j})u_{m,j}}{\Delta x^{2}}-\frac{u_{m,j-1}-2u_{m,j}+u_{m,j+1}}{\Delta y^{2}}=f_{m,j}+\frac{2}{\Delta x}\varphi_{2}(y_{j}),1\leqslant j\leqslant n-1,\\ -\frac{u_{i-1,0}-2u_{i,0}+u_{i+1,0}}{\Delta x^{2}}-\frac{2u_{i,1}-2(1+\Delta y\lambda_{i,0})u_{i,0}}{\Delta y^{2}}=f_{i,0}-\frac{2}{\Delta y}\Psi_{1}(x_{i}),1\leqslant i\leqslant m-1,\\ -\frac{u_{i-1,n}-2u_{i,n}+u_{i+1,n}}{\Delta x^{2}}-\frac{2u_{i,n-1}-2(1+\Delta y\lambda_{i,n})u_{i,n}}{\Delta y^{2}}=f_{i,n}+\frac{2}{\Delta y}\Psi_{2}(x_{i}),1\leqslant i\leqslant m-1 \end{matrix}\right.

        至此我们共有(m+1)x(n+1)个待求量u_{i,j},0\leqslant i\leqslant m,0\leqslant j\leqslant n,而现有(m-1)x(n-1)个关于内节点的方程,2(n-1)个关于左、右边界上的节点(不含端点)的方程及2(m-1)个关于上、下边界上的节点(也不含端点)的方程,还需要补充

(m+1)(n+1)-(m-1)(n-1)-2(n-1)-2(m-1)=4

个方程,也就是关于矩形区域的4个顶点的方程。为此,设公式(2)中第1式对i=0,j=0成立,即

-\frac{u_{-1,0}-2u_{0,0}+u_{1,0}}{\Delta x^{2}}-\frac{u_{0,-1}-2u_{0,0}+u_{0,1}}{\Delta y^{2}}=f_{0,0} \; (3)

        然后再从公式(2)的第2式和第4式中分别解出u_{-1,0}=u_{1,0}-2\Delta x(\varphi_{1}(y_{0})+\lambda_{0,0}u_{0,0})u_{0,-1}=u_{0,1}-2\Delta y(\Psi_{1}(x_{0})+\lambda_{0,0}u_{0,0}),代入公式(3)中,得到\Omega左下顶点处的方程:

-\frac{2u_{1,0}-2(1+\Delta x\lambda_{0,0})u_{0,0}}{\Delta x^{2}}-\frac{2u_{0,1}-2(1+\Delta y\lambda_{0,0})u_{0,0}}{\Delta y^{2}}=f_{0,0}-\frac{2}{\Delta x}\varphi_{1}(y_{0})-\frac{2}{\Delta y}\Psi_{1}(x_{0})

上式可以简化为

-\frac{2u_{1,0}}{\Delta x^{2}}+[\frac{2(1+\Delta x\lambda_{0,0})}{\Delta x^{2}}+\frac{2(1+\Delta y\lambda_{0,0})}{\Delta y^{2}}]u_{0,0}-\frac{2u_{0,1}}{\Delta y^{2}}=f_{0,0}-\frac{2}{\Delta x}\varphi_{1}(y_{0})-\frac{2}{\Delta y}\Psi_{1}(x_{0})

        同样,设公式(2)中第1式分别对i=m,j=0i=0,j=ni=m,j=n成立,然后再从公式(2)中第3式和第4式解出u_{m+1,0},u_{m,-1},u_{-1,n},u_{0,n+1},u_{m+1,n},u_{m,n+1}分别代入刚才得到的3个方程,就得到\Omega的右下顶点、左上顶点、和右上顶点处的方程,这样,我们就有了完整的处理带导数边界条件的椭圆型方程的数值格式:

\left\{\begin{matrix} -\frac{u_{i,j-1}}{\Delta y^{2}}-\frac{u_{i-1,j}}{\Delta x^{2}}+2[\frac{1}{\Delta x^{2}}+\frac{1}{\Delta y^{2}}]u_{i,j}-\frac{u_{i+1,j}}{\Delta x^{2}}-\frac{u_{i,j+1}}{\Delta y^{2}}=f_{i,j},1\leqslant i\leqslant m-1,1\leqslant j\leqslant n-1,\\ -\frac{u_{0,j-1}}{\Delta y^{2}}+[\frac{2(1+\Delta x\lambda_{0,j})}{\Delta x^{2}}+\frac{2}{\Delta y^{2}}]u_{0,j}-\frac{2u_{1,j}}{\Delta x^{2}}-\frac{u_{0,j+1}}{\Delta y^{2}}=f_{0,j}-\frac{2}{\Delta x}\varphi_{1}(y_{j}),1\leqslant j\leqslant n-1,\\ -\frac{u_{m,j-1}}{\Delta y^{2}}-\frac{2u_{m-1,j}}{\Delta x^{2}}+[\frac{2(1+\Delta x\lambda_{m,j})}{\Delta x^{2}}+\frac{2}{\Delta y^{2}}]u_{m,j}-\frac{u_{m,j+1}}{\Delta y^{2}}=f_{m,j}+\frac{2}{\Delta x}\varphi_{2}(y_{j}),1\leqslant j\leqslant n-1,\\ -\frac{u_{i-1,0}}{\Delta x^{2}}+[\frac{2}{\Delta x^{2}}+\frac{2(1+\Delta y\lambda_{i,0})}{\Delta y^{2}}]u_{i,0}-\frac{u_{i+1,0}}{\Delta x^{2}}-\frac{2u_{i,1}}{\Delta y^{2}}=f_{i,0}-\frac{2}{\Delta y}\Psi_{1}(x_{i}),1\leqslant i\leqslant m-1,\\ -\frac{u_{i-1,n}}{\Delta x^{2}}-\frac{2u_{i,n-1}}{\Delta y^{2}}+[\frac{2}{\Delta x^{2}}+\frac{2(1+\Delta y\lambda_{i,n})}{\Delta y^{2}}]u_{i,n-\frac{u_{i+1,n}}{\Delta x^{2}}}=f_{i,n}+\frac{2}{\Delta y}\Psi_{2}(x_{i}),1\leqslant i\leqslant m-1,\\ [\frac{2(1+\Delta x\lambda_{0,0})}{\Delta x^{2}}+\frac{2(1+\Delta y \lambda_{0,0})}{\Delta y^{2}}]u_{0,0}-\frac{2u_{1,0}}{\Delta x^{2}}-\frac{2u_{0,1}}{\Delta y^{2}}=f_{0,0}-\frac{2}{\Delta x}\varphi_{1}(y_{0})-\frac{2}{\Delta y}\Psi_{1}(x_{0}),\\ -\frac{2u_{m-1,0}}{\Delta x^{2}}+[\frac{2(1+\Delta x\lambda_{m,0})}{\Delta x^{2}}+\frac{2(1+\Delta y \lambda_{m,0})}{\Delta y^{2}}]u_{m,0}-\frac{2u_{m,1}}{\Delta y^{2}}=f_{m,0}+\frac{2}{\Delta x}\varphi_{2}(y_{0})-\frac{2}{\Delta y}\Psi_{1}(x_{m}),\\ -\frac{2u_{0,n-1}}{\Delta y^{2}}+[\frac{2(1+\Delta x\lambda_{0,n})}{\Delta x^{2}}+\frac{2(1+\Delta y\lambda_{0,n})}{\Delta y^{2}}]u_{0,n}-\frac{2u_{1,n}}{\Delta x^{2}}=f_{0,n}-\frac{2}{\Delta x}\varphi_{1}(y_{n})+\frac{2}{\Delta y}\Psi_{2}(x_{0}),\\ -\frac{2u_{m,n-1}}{\Delta y^{2}}-\frac{2u_{m-1,n}}{\Delta x^{2}}+[\frac{2(1+\Delta x\lambda_{m,n})}{\Delta x^{2}}+\frac{2(1+\Delta y\lambda_{m,n})}{\Delta y^{2}}]u_{m,n}=f_{m,n}+\frac{2}{\Delta x}\varphi_{2}(y_{n})+\frac{2}{\Delta y}\Psi_{2}(x_{m}). \end{matrix}\right.

        记\beta =\frac{1}{\Delta x^{2}},\gamma=\frac{1}{\Delta y^{2}},\alpha=2(\frac{1}{\Delta x^{2}}+\frac{1}{\Delta y^{2}}),\xi=\frac{2}{\Delta x},\eta=\frac{2}{\Delta y},有

\begin{Bmatrix} -\gamma u_{i,j-1}-\beta u_{i-1,j}+\alpha u_{i,j}-\beta u_{i+1,j}-\gamma u_{i,j+1}=f_{i,j},1\leqslant i\leqslant m-1,1\leqslant j\leqslant n-1,\\ -\gamma u_{0,j-1}+[\alpha+\xi \lambda_{0,j}]u_{0,j}-2\beta u_{1,j}-\gamma u_{0,j+1}=f_{0,j}-\xi\varphi_{1}(y_{j}),1\leqslant j\leqslant n-1,\\ -\gamma u_{m,j-1}-2\beta u_{m-1,j}+[\alpha+\xi\lambda_{m,j}]u_{m,j}-\gamma u_{m,j+1}=f_{m,j}+\xi \varphi_{2}(y_{j}),1\leqslant j\leqslant n-1,\\ -\beta u_{i-1,0}+[\alpha+\eta\lambda_{i,0}]u_{i,0}-\beta u_{i+1,0}-2\gamma u_{i,1}=f_{i,0}-\eta\Psi_{1}(x_{i}),1\leqslant i\leqslant m-1,\\ -2\gamma u_{i,n-1}-\beta u_{i-1,n}+[\alpha+\eta\lambda_{i,n}]u_{i,n}-\beta u_{i+1,n}=f_{i,n}+\eta\Psi_{2}(x_{i}),1\leqslant i\leqslant m-1,\\ [\alpha+(\xi+\eta)\lambda_{0,0}]u_{0,0}-2\beta u_{1,0}-2\gamma u_{0,1}=f_{0,0}-\xi\varphi_{1}(y_{0})-\eta\Psi_{1}(x_{0}),\\ -2\beta u_{m-1,0}+[\alpha+(\xi+\eta)\lambda_{m,0}]u_{m,0}-2\gamma u_{m,1}=f_{m,0}+\xi\varphi_{2}(y_{0})-\eta\Psi_{1}(x_{m}),\\ -2\gamma u_{0,n-1}+[\alpha+(\xi+\eta)\lambda_{0,n}]u_{0,n}-2\beta u_{1,n}=f_{0,n}-\xi \varphi_{1}(y_{n})+\eta\Psi_{2}(x_{0}),\\ -2\gamma u_{m,n-1}-2\beta u_{m-1,n}+[\alpha+(\xi+\eta)\lambda_{m,n}]u_{m,n}=f_{m,n}+\xi\varphi_{2}(y_{n})+\eta\Psi_{2}(x_{m}). \end{Bmatrix}\;(4)

        为计算分别,可将上述方程组写成矩阵形式。

        首先,将公式(4)中第4、6、7式写成:

\begin{pmatrix} \alpha+(\xi+\eta)\lambda_{0,0} & -2\beta & & & & & \\ -\beta & \alpha+\eta\lambda_{1,0} & -\beta & & & & \\ & -\beta & \alpha+\eta\lambda_{2,0} & -\beta & & & \\ & & \ddots & \ddots & \ddots & & \\ & & & -\beta & \alpha+\eta\lambda_{m-2,0} & -\beta & \\ & & & & -\beta & \alpha+\eta\lambda_{m-1,0} & -\beta\\ & & & & & -2\beta & \alpha+(\xi+\eta)\lambda_{m,0} \end{pmatrix}+\begin{pmatrix} u_{0,0}\\ u_{1,0}\\ u_{2,0}\\ \vdots\\ u_{m-2,0}\\ u_{m-1,0}\\ u_{m,0} \end{pmatrix}+\begin{pmatrix} -2\gamma & & & & & & \\ & -2\gamma & & & & & \\ & & -2\gamma & & & & \\ & & & -2\gamma & & & \\ & & & & -2\gamma & & \\ & & & & & -2\gamma & \\ & & & & & & -2\gamma \end{pmatrix}\begin{pmatrix} u_{0,1}\\ u_{1,1}\\ u_{2,1}\\ \vdots\\ u_{m-2,1}\\ u_{m-1,1}\\ u_{m,1} \end{pmatrix}=\begin{pmatrix} f_{0,0}-\eta\Psi_{1}(x_{0})-\xi\varphi_{1}(y_{0})\\ f_{1,0}-\eta\Psi_{1}(x_{1})\\ f_{2,0}-\eta\Psi_{1}(x_{2})\\ \vdots\\ f_{m-2,0}-\eta\Psi_{1}(x_{m-2})\\ f_{m-1,0}-\eta\Psi_{1}(x_{m-1})\\ f_{m,0}-\eta\Psi_{1}(x_{m})-\xi\varphi_{2}(y_{0}) \end{pmatrix} \; (5)

         上面的式子可以简记为

C\mathbf{u_{0}}+2A\mathbf{u}_{1}=\mathbf{f}_{0}

其中,A=-\gamma II为m+1阶单位矩阵,且C为公式(5)最左端的三对角矩阵,\mathbf{f}_{0}为公式(5)右端的向量,\mathbf{u}_{j}=(u_{0,j},u_{1,j},\cdot\cdot\cdot ,u_{m-1,j},u_{m,j})^{T},0\leqslant j\leqslant n。接着,公式(4)中的第1,2,3式可以写成

\begin{pmatrix} -\gamma & & & & & & \\ & -\gamma & & & & & \\ & & -\gamma & & & & \\ & &\ddots &\ddots &\ddots & & \\ & & & & -\gamma & & \\ & & & & & -\gamma & \\ & & & & & & -\gamma \end{pmatrix}\begin{pmatrix} u_{0,j-1}\\ u_{1,j-1}\\ u_{2,j-1}\\ \vdots\\ u_{m-2,j-1}\\ u_{m-1,j-1}\\ u_{m,j-1} \end{pmatrix}+\begin{pmatrix} \alpha+\xi\lambda_{0,j} & -2\beta & & & & & \\ -\beta & \alpha & -\beta & & & & \\ & -\beta & \alpha & -\beta & & & \\ & & \ddots & \ddots & \ddots & & \\ & & & -\beta & \alpha & -\beta & \\ & & & & -\beta & \alpha & -\beta\\ & & & & & -2\beta & \alpha+\xi\lambda_{m,j} \end{pmatrix}\begin{pmatrix} u_{0,j}\\ u_{1,j}\\ u_{2,j}\\ \vdots\\ u_{m-2,j}\\ u_{m-1,j}\\ u_{m,j} \end{pmatrix}+\begin{pmatrix} -\gamma & & & & & & \\ & -\gamma & & & & & \\ & & -\gamma & & & & \\ & &\ddots &\ddots &\ddots & & \\ & & & & -\gamma & & \\ & & & & & -\gamma & \\ & & & & & & -\gamma \end{pmatrix}\begin{pmatrix} u_{0,j+1}\\ u_{1,j+1}\\ u_{2,j+1}\\ \vdots\\ u_{m-2,j+1}\\ u_{m-1,j+1}\\ u_{m,j+1} \end{pmatrix}=\begin{pmatrix} f_{0,j}-\xi \varphi_{1}(y_{j})\\ f_{1,j}\\ f_{2,j}\\ \vdots\\ f_{m-2,j}\\ f_{m-1,j}\\ f_{m,j}+\xi \varphi_{2}(y_{j}) \end{pmatrix}\;(6),1\leqslant j\leqslant n-1

         该式可以简记为

A\mathbf{u}_{j-1}+B_{j}\mathbf{u}_{j}+A\mathbf{u}_{j+1}=\mathbf{f}_{j},1\leqslant j\leqslant n-1

其中,B_{j}为公式(6)中的三对角矩阵,\mathbf{f}_{j}为公式(6)右端的向量。最后,公式(4)的第5、8、9式可以写成

\begin{pmatrix} -2\gamma & & & & & & \\ & -2\gamma & & & & & \\ & & -2\gamma & & & & \\ & & \ddots& \ddots& \ddots& & \\ & & & & -2\gamma & & \\ & & & & & -2\gamma & \\ & & & & & & -2\gamma \end{pmatrix}\begin{pmatrix} u_{0,n-1}\\ u_{1,n-1}\\ u_{2,n-1}\\ \vdots\\ u_{m-2,n-1}\\ u_{m-1,n-1}\\ u_{m,n-1} \end{pmatrix}+\begin{pmatrix} \alpha+(\xi+\eta)\lambda_{0,n} & -2\beta & & & & & \\ -\beta & \alpha+\eta\lambda_{1,n} & -\beta & & & & \\ & -\beta & \alpha+\eta\lambda_{2,n} & -\beta & & & \\ & & \ddots & \ddots & \ddots & & \\ & & & -\beta & \alpha+\eta\lambda_{m-2,n} & -\beta & \\ & & & & -\beta & \alpha+\eta\lambda_{m-1,n} & -\beta\\ & & & & & -2\beta & \alpha+(\xi+\eta)\lambda_{m,n} \end{pmatrix}\begin{pmatrix} u_{0,n}\\ u_{1,n}\\ u_{2,n}\\ \vdots\\ u_{m-2,n}\\ u_{m-1,n}\\ u_{m,n} \end{pmatrix}==\begin{pmatrix} f_{0,n}+\eta\Psi_{1}(x_{0})-\xi\varphi_{1}(y_{n})\\ f_{1,n}+\eta\Psi_{2}(x_{1})\\ f_{2,n}-\eta\Psi_{2}(x_{2})\\ \vdots\\ f_{m-2,n}+\eta\Psi_{2}(x_{m-2})\\ f_{m-1,n}+\eta\Psi_{2}(x_{m-1})\\ f_{m,n}+\eta\Psi_{2}(x_{m})+\xi\varphi_{2}(y_{n}) \end{pmatrix} \; (7)

        该式可以简记为

2A\mathbf{u}_{n-1}+D\mathbf{u}_{n}=\mathbf{f}_{n}

其中,D为公式(7)中的三对角矩阵,\mathbf{f}_{n}为公式(7)右端的向量。于是,由公式(5)、(6)、(7)可得到公式(4)写成块三对角矩阵为

\begin{pmatrix} C & 2A & & & & & \\ A & B_{1} & A & & & & \\ & A & B_{2} & A & & & \\ & & \ddots & \ddots & \ddots & & \\ & & & A & B_{m-2} & A & \\ & & & & A & B_{m-1} & A\\ & & & & & 2A & D \end{pmatrix}\begin{pmatrix} \mathbf{u}_{0}\\ \mathbf{u}_{1}\\ \mathbf{u}_{2}\\ \vdots\\ \mathbf{u}_{m-2}\\ \mathbf{u}_{m-1}\\ \mathbf{u}_{m} \end{pmatrix}=\begin{pmatrix} \mathbf{f}_{0}\\ \mathbf{f}_{1}\\ \mathbf{f}_{2}\\ \vdots\\ \mathbf{f}_{m-2}\\ \mathbf{f}_{m-1}\\ \mathbf{f}_{m} \end{pmatrix}\;(8)

        采用Gauss-Seidel迭代法求解公式(8)。

三、算例实现

        用差分格式-公式(8)求解椭圆型方程混合边值问题:

\left\{\begin{matrix} -(\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}})=(\pi^{2}-1)e^{x}sin(\pi y),0<x<2,0<y<1,\\ (\frac{\partial u(x,y)}{\partial x}-(x^{2}+y^{2})u)|_{(0,y)}=(1-y^{2})sin(\pi y),0\leqslant y\leqslant 1,\\ (\frac{\partial u(x,y)}{\partial x}+(x^{2}+y^{2})u)|_{(2,y)}=(5+y^{2})e^{2}sin(\pi y),0\leqslant y\leqslant 1,\\ (\frac{\partial u(x,y)}{\partial y}-(x^{2}+y^{2})u)|_{(x,c)}=\pi e^{x},0\leqslant x\leqslant 1,\\ (\frac{\partial u(x,y)}{\partial y}+(x^{2}+y^{2})u)|_{(x,d)}=-\pi e^{x},0\leqslant x\leqslant 1 \end{matrix}\right.

已知精确解为u(x,y)=e^{x}sin(\pi y)。分别取步长为\Delta x=\Delta y=\frac{1}{16}\Delta x=\Delta y=\frac{1}{32},输出6个节点(0.5i,0.25)(0.5i,0.50),i=1,2,3处的数值解,并给出误差,要求在各节点处最大误差的迭代误差限为0.5\times10^{-10}

代码如下:


#include <cmath>
#include <stdlib.h>
#include <stdio.h>
#define pi 3.14159265359int main(int argc, char*argv[])
{int m,n,i,j,k;double xa,xb,ya,yb,dx,dy,alpha,beta,gamma,maxerr;double *x,*y,**u,**v,**lambda,*d,kexi,eta,temp;double f(double x, double y);double lambda_function(double x, double y);double phi1(double y);double phi2(double y);double psi1(double x);double psi2(double x);double exact(double x, double y);xa=0.0;xb=2.0;ya=0.0;yb=1.0;m=64;n=32;printf("m=%d,n=%d.\n",m,n);dx=(xb-xa)/m;dy=(yb-ya)/n;beta=1.0/(dx*dx);gamma=1.0/(dy*dy);alpha=2*(beta+gamma);kexi=2.0/dx;eta=2.0/dy;x=(double*)malloc(sizeof(double)*(m+1));for(i=0;i<=m;i++)x[i]=xa+i*dx;y=(double*)malloc(sizeof(double)*(n+1));for(j=0;j<=n;j++)y[j]=ya+j*dy;u=(double**)malloc(sizeof(double*)*(m+1));v=(double**)malloc(sizeof(double*)*(m+1));lambda=(double**)malloc(sizeof(double*)*(m+1));for(i=0;i<=m;i++){u[i]=(double*)malloc(sizeof(double)*(n+1));v[i]=(double*)malloc(sizeof(double)*(n+1));lambda[i]=(double*)malloc(sizeof(double)*(n+1));}for(i=0;i<=m;i++){for(j=0;j<=n;j++){u[i][j]=0.0;v[i][j]=0.0;lambda[i][j]=lambda_function(x[i],y[j]);}}d=(double*)malloc(sizeof(double)*(m+1));k=0;do{maxerr=0.0;for(i=0;i<=m;i++)d[i]=f(x[i],y[0])-eta*psi1(x[i]);d[0]=d[0]-kexi*phi1(y[0]);d[m]=d[m]+kexi*phi2(y[0]);v[0][0]=(d[0]+2*gamma*u[0][1]+2*beta*u[1][0])/(alpha+(kexi+eta)*lambda[0][0]);for(i=1;i<m;i++)v[i][0]=(d[i]+2*gamma*u[i][1]+beta*(v[i-1][0]+u[i+1][0]))/(alpha+eta*lambda[i][0]);v[m][0]=(d[m]+2*gamma*u[m][1]+2*beta*v[m-1][0])/(alpha+(kexi+eta)*lambda[m][0]);for(j=1;j<n;j++){for(i=0;i<=m;i++){d[i]=f(x[i],y[j]);}d[0]=d[0]-kexi*phi1(y[j]);d[m]=d[m]+kexi*phi2(y[j]);v[0][j]=(d[0]+gamma*(u[0][j+1]+v[0][j-1])+2*beta*u[1][j])/(alpha+kexi*lambda[0][j]);for(i=1;i<m;i++){v[i][j]=(d[i]+gamma*(v[i][j-1]+u[i][j+1])+beta*(v[i-1][j]+u[i+1][j]))/alpha;}v[m][j]=(d[m]+gamma*(v[m][j-1]+u[m][j+1])+2*beta*v[m-1][j])/(alpha+kexi*lambda[m][j]);}for(i=0;i<=m;i++)d[i]=f(x[i],y[n])+eta*psi2(x[i]);d[0]=d[0]-kexi*phi1(y[n]);d[m]=d[m]+kexi*phi2(y[n]);v[0][n]=(d[0]+2*gamma*v[0][n-1]+2*beta*u[1][n])/(alpha+(kexi+eta)*lambda[0][n]);for(i=1;i<m;i++)v[i][n]=(d[i]+2*gamma*v[i][n-1]+beta*(v[i-1][n]+u[i+1][n]))/(alpha+eta*lambda[i][n]);v[m][n]=(d[m]+2*gamma*v[m][n-1]+2*beta*v[m-1][n])/(alpha+(kexi+eta)*lambda[m][n]);for(i=0;i<=m;i++){for(j=0;j<=n;j++){temp=fabs(u[i][j]-v[i][j]);if(temp>maxerr)maxerr=temp;u[i][j]=v[i][j];}}k=k+1;}while((maxerr>0.5*1e-10)&&(k<=1e+8));printf("k=%d\n",k);k=m/4;for(i=k;i<m;i=i+k){printf("(%.2f,0.25), y=%f, err=%.4e.\n",x[i],u[i][n/4],fabs(exact(x[i],y[n/4])-u[i][n/4]));}k=m/4;for(i=k;i<m;i=i+k){printf("(%.2f,0.50), y=%f, err=%.4e.\n",x[i],u[i][n/2],fabs(exact(x[i],y[n/2])-u[i][n/2]));}for(i=0;i<=m;i++){free(u[i]);free(v[i]);free(lambda[i]);}free(u);free(v);free(lambda);free(x);free(y);free(d);return 0;
}double f(double x, double y)
{return (pi*pi-1)*exp(x)*sin(pi*y);
}
double lambda_function(double x, double y)
{return x*x+y*y;
}
double phi1(double y)
{return (1-y*y)*sin(pi*y);
}
double phi2(double y)
{return (5.0+y*y)*exp(2.0)*sin(pi*y);
}
double psi1(double x)
{return pi*exp(x);
}
double psi2(double x)
{return -pi*exp(x);
}
double exact(double x, double y)
{return exp(x)*sin(pi*y);
}

\Delta x=\Delta y=\frac{1}{16}时,计算结果如下:

m=32,n=16.
k=4323
(0.50,0.25), y=1.152179, err=1.3643e-02.
(1.00,0.25), y=1.911016, err=1.1100e-02.
(1.50,0.25), y=3.162159, err=6.8738e-03.
(0.50,0.50), y=1.638607, err=1.0115e-02.
(1.00,0.50), y=2.711255, err=7.0265e-03.
(1.50,0.50), y=4.479936, err=1.7526e-03.

\Delta x=\Delta y=\frac{1}{32}时,计算结果如下:

m=64,n=32.
k=16007
(0.50,0.25), y=1.162412, err=3.4097e-03.
(1.00,0.25), y=1.919341, err=2.7745e-03.
(1.50,0.25), y=3.167313, err=1.7193e-03.
(0.50,0.50), y=1.646193, err=2.5286e-03.
(1.00,0.50), y=2.716524, err=1.7575e-03.
(1.50,0.50), y=4.481249, err=4.3972e-04.

四、结论

        从计算结果可知,当步长减小为1/2时,误差减小为原来的1/4,可见该数值格式是二阶收敛的。

相关文章:

偏微分方程算法之混合边界条件下的差分法

目录 一、研究目标 二、理论推导 三、算例实现 四、结论 一、研究目标 我们在前几节中介绍了Poisson方程的边值问题&#xff0c;接下来对椭圆型偏微分方程的混合边值问题进行探讨&#xff0c;研究对象为&#xff1a; 其中&#xff0c;为矩形区域&#xff0c;为上的连续函数…...

apollo资料整理

Application X: Application X Apollo: Apollo 自动驾驶开放平台 Cyber RT API tutorial — Cyber RT Documents documentation Cyber RT API tutorial — Cyber RT Documents documentation GitHub - daohu527/dig-into-apollo: Apollo notes (Apollo学习笔记) - Apollo l…...

森林消防新利器:高扬程水泵的革新与应用/恒峰智慧科技

随着全球气候变化的加剧&#xff0c;森林火灾的频发已成为威胁生态安全的重要问题。在森林消防工作中&#xff0c;高效、快速的水源供给设备显得尤为重要。近年来&#xff0c;高扬程水泵的广泛应用&#xff0c;为森林消防工作带来了新的希望与突破。 一、高扬程水泵的技术优势 …...

Microsoft Universal Print 与 SAP 集成教程

引言 从 SAP 环境打印是许多客户的要求。例如数据列表打印、批量打印或标签打印。此类生产和批量打印方案通常使用专用硬件、驱动程序和打印解决方案来解决。 Microsoft Universal Print 是一种基于云的打印解决方案&#xff0c;它允许组织以集中化的方式管理打印机和打印机驱…...

VBA在Excel中字母、数字的相互转化

VBA在Excel中字母、数字的相互转化 字母转数字的方法 数字转字母的方法 众所周知,Excel表中的行以数字展示,列用字母展示,如下图: 编程时,很多时候需要将列的字母转变为数字使用,如cells(num1,num2).value等,不知大家是怎么将字母转化为数字的,Excel是否有其他方式…...

【C语言】——联合体与枚举

【C语言】——联合体与枚举 一、联合体1.1、联合体类型的声明1.2、联合体的特点1.3、相同成员的结构体和联合体对比1.4、联合体的大小计算1.5、联合体的应用举例 二、枚举2.1、枚举类型的声明2.2、枚举类型的优点 一、联合体 1.1、联合体类型的声明 联合体也叫做共用体   与…...

java线上问题排查之内存分析(三)

java线上问题排查之内存分析 使用top命令 top命令显示的结果列表中&#xff0c;会看到%MEM这一列&#xff0c;这里可以看到你的进程可能对内存的使用率特别高。以查看正在运行的进程和系统负载信息&#xff0c;包括cpu负载、内存使用、各个进程所占系统资源等。 2.用jstat命令…...

中电金信:金Gien乐道 | 4月要闻速览,精彩再回顾

中国电子党组副书记、总经理李立功一行调研中电金信 4月10日&#xff0c;中国电子党组副书记、总经理李立功一行赴中电金信进行调研&#xff0c;深入听取了中电金信经营发展情况、研发工作及“源启”行业数字底座平台的汇报&#xff0c;并参观了公司展厅和科技研发场所&#xf…...

Java将文件目录转成树结构

在实际开发中经常会遇到返回树形结构的场景&#xff0c;特别是在处理文件系统或者是文件管理系统中。下面就介绍一下怎么将文件路径转成需要的树形结构。 在Java中&#xff0c;将List<String>转换成树状结构&#xff0c;需要定义一个树节点类&#xff08;TreeNode&#…...

硬件工程师必读:10条职业发展黄金法则!

在快速发展的科技时代&#xff0c;硬件工程师作为推动技术创新和产业升级的重要力量&#xff0c;其职业发展之路既充满挑战也蕴含无限机遇。为了在这条道路上稳步前行&#xff0c;我们首先需要了解硬件产品的研发流程。 在这个过程中&#xff0c;公司内的每个岗位都发挥着不可或…...

Redis是什么? 日常运维 Redis 需要注意什么 ? 怎么降低Redis 内存使用 节省内存?

你的项目或许已经使用 Redis 很长时间了&#xff0c;但在使用过程中&#xff0c;你可能还会或多或少地遇到以下问题&#xff1a; 我的 Redis 内存为什么增长这么快&#xff1f;为什么我的 Redis 操作延迟变大了&#xff1f;如何降低 Redis 故障发生的频率&#xff1f;日常运维…...

【Android项目】“追茶到底”项目介绍

没有多的介绍&#xff0c;这里只是展示我的项目效果&#xff0c;后面会给出具体的代码实现。 一、用户模块 1、注册&#xff08;第一次登陆的话需要先注册账号&#xff09; 2、登陆&#xff08;具有记住最近登录用户功能&#xff09; 二、点单模块 1、展示饮品列表 2、双向联动…...

机试:进制转换问题

十进制转任意进制 简单回忆一下十进制我们是怎么转换成二进制的&#xff08;短除法&#xff09;&#xff1a; 我们会将十进制数不断的进行除2操作&#xff0c;并且记录下每一次的余数&#xff08;这个余数就是我们最终求的二进制数的组成部分&#xff09;。 以下以12D举例&a…...

目标检测实战(十五): 使用YOLOv7完成对图像的目标检测任务(从数据准备到训练测试部署的完整流程)

文章目录 一、目标检测介绍二、YOLOv7介绍三、源码/论文获取四、环境搭建4.1 环境检测 五、数据集准备六、 模型训练七、模型验证八、模型测试九、错误总结9.1 错误1-numpy jas mp attribute int9.2 错误2-测试代码未能跑出检测框9.3 错误3- Command git tag returned non-zero…...

github中fasttext库README官文文档翻译

参考链接&#xff1a;fastText/python/README.md at main facebookresearch/fastText (github.com) fastText模块介绍 fastText 是一个用于高效学习单词表述和句子分类的库。在本文档中&#xff0c;我们将介绍如何在 python 中使用 fastText。 环境要求 fastText 可在现代 …...

WouoUIPagePC端实现

WouoUIPagePC端实现 WouoUIPage是一个与硬件平台无关&#xff0c;纯C语言的UI库&#xff08;目前只能应用于128*64的单色OLED屏幕上&#xff0c;后期会改进&#xff0c;支持更多尺寸&#xff09;。因此&#xff0c;我们可以在PC上实现它&#xff0c;本文就以在PC上使用 VScode…...

W801学习笔记十九:古诗学习应用——下

经过前两章的内容&#xff0c;背唐诗的功能基本可以使用了。然而&#xff0c;仅有一种模式未免显得过于单一。因此&#xff0c;在本章中对其进行扩展&#xff0c;增加几种不同的玩法&#xff0c;并且这几种玩法将采用完全不同的判断方式。 玩法一&#xff1a;三分钟限时挑战—…...

类加载器ClassLoad-jdk1.8

类加载器ClassLoad-jdk1.8 1. 类加载器的作用2. 类加载器的种类&#xff08;JDK8&#xff09;3. jvm内置类加载器如何搜索加载类--双亲委派模型4. 如何打破双亲委派模型--自定义类加载器5. 自定义一个类加载器5.1 为什么需要自定义类加载器5.2 自定义一个类加载器 6. java代码加…...

24年最新AI数字人简单混剪

24年最新AI数字人简单混剪 网盘自动获取 链接&#xff1a;https://pan.baidu.com/s/1lpzKPim76qettahxvxtjaQ?pwd0b8x 提取码&#xff1a;0b8x...

免备案香港主机会影响网站收录?

免备案香港主机会影响网站收录?前几天遇到一个做电子商务的朋友说到这个使用免备案香港主机的完整会不会影响网站的收录问题&#xff0c;这个问题也是站长关注较多的问题之一。小编查阅了百度官方规则说明&#xff0c;应该属于比较全面的。下面小编给大家介绍一下使用免备案香…...

Ubuntu系统下交叉编译openssl

一、参考资料 OpenSSL&&libcurl库的交叉编译 - hesetone - 博客园 二、准备工作 1. 编译环境 宿主机&#xff1a;Ubuntu 20.04.6 LTSHost&#xff1a;ARM32位交叉编译器&#xff1a;arm-linux-gnueabihf-gcc-11.1.0 2. 设置交叉编译工具链 在交叉编译之前&#x…...

在鸿蒙HarmonyOS 5中实现抖音风格的点赞功能

下面我将详细介绍如何使用HarmonyOS SDK在HarmonyOS 5中实现类似抖音的点赞功能&#xff0c;包括动画效果、数据同步和交互优化。 1. 基础点赞功能实现 1.1 创建数据模型 // VideoModel.ets export class VideoModel {id: string "";title: string ""…...

安宝特方案丨XRSOP人员作业标准化管理平台:AR智慧点检验收套件

在选煤厂、化工厂、钢铁厂等过程生产型企业&#xff0c;其生产设备的运行效率和非计划停机对工业制造效益有较大影响。 随着企业自动化和智能化建设的推进&#xff0c;需提前预防假检、错检、漏检&#xff0c;推动智慧生产运维系统数据的流动和现场赋能应用。同时&#xff0c;…...

JVM垃圾回收机制全解析

Java虚拟机&#xff08;JVM&#xff09;中的垃圾收集器&#xff08;Garbage Collector&#xff0c;简称GC&#xff09;是用于自动管理内存的机制。它负责识别和清除不再被程序使用的对象&#xff0c;从而释放内存空间&#xff0c;避免内存泄漏和内存溢出等问题。垃圾收集器在Ja…...

第一篇:Agent2Agent (A2A) 协议——协作式人工智能的黎明

AI 领域的快速发展正在催生一个新时代&#xff0c;智能代理&#xff08;agents&#xff09;不再是孤立的个体&#xff0c;而是能够像一个数字团队一样协作。然而&#xff0c;当前 AI 生态系统的碎片化阻碍了这一愿景的实现&#xff0c;导致了“AI 巴别塔问题”——不同代理之间…...

听写流程自动化实践,轻量级教育辅助

随着智能教育工具的发展&#xff0c;越来越多的传统学习方式正在被数字化、自动化所优化。听写作为语文、英语等学科中重要的基础训练形式&#xff0c;也迎来了更高效的解决方案。 这是一款轻量但功能强大的听写辅助工具。它是基于本地词库与可选在线语音引擎构建&#xff0c;…...

短视频矩阵系统文案创作功能开发实践,定制化开发

在短视频行业迅猛发展的当下&#xff0c;企业和个人创作者为了扩大影响力、提升传播效果&#xff0c;纷纷采用短视频矩阵运营策略&#xff0c;同时管理多个平台、多个账号的内容发布。然而&#xff0c;频繁的文案创作需求让运营者疲于应对&#xff0c;如何高效产出高质量文案成…...

使用LangGraph和LangSmith构建多智能体人工智能系统

现在&#xff0c;通过组合几个较小的子智能体来创建一个强大的人工智能智能体正成为一种趋势。但这也带来了一些挑战&#xff0c;比如减少幻觉、管理对话流程、在测试期间留意智能体的工作方式、允许人工介入以及评估其性能。你需要进行大量的反复试验。 在这篇博客〔原作者&a…...

何谓AI编程【02】AI编程官网以优雅草星云智控为例建设实践-完善顶部-建立各项子页-调整排版-优雅草卓伊凡

何谓AI编程【02】AI编程官网以优雅草星云智控为例建设实践-完善顶部-建立各项子页-调整排版-优雅草卓伊凡 背景 我们以建设星云智控官网来做AI编程实践&#xff0c;很多人以为AI已经强大到不需要程序员了&#xff0c;其实不是&#xff0c;AI更加需要程序员&#xff0c;普通人…...

GB/T 43887-2024 核级柔性石墨板材检测

核级柔性石墨板材是指以可膨胀石墨为原料、未经改性和增强、用于核工业的核级柔性石墨板材。 GB/T 43887-2024核级柔性石墨板材检测检测指标&#xff1a; 测试项目 测试标准 外观 GB/T 43887 尺寸偏差 GB/T 43887 化学成分 GB/T 43887 密度偏差 GB/T 43887 拉伸强度…...