偏微分方程算法之二阶双曲型方程交替方向隐格式

2024-04-25 02:36

本文主要是介绍偏微分方程算法之二阶双曲型方程交替方向隐格式,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

目录

一、研究目标

二、理论推导

2.1 显差分格式

2.2 交替方向隐格式

三、算例实现

四、结论


一、研究目标

        本节的研究对象为如下的双曲型偏微分方程初边值问题:

\left\{\begin{matrix} \frac{\partial^{2}u(x,y,t)}{\partial t^{2}}-(\frac{\partial^{2}u(x,y,t)}{\partial x^{2}}+\frac{\partial^{2}u(x,y,t)}{\partial y^{2}})=f(x,y,t),(x,y)\in \Omega=[0,a]\times [0,b],0<t\leqslant T,\\ u(x,y,0)=\varphi(x,y),\frac{\partial u}{\partial t}(x,y,0)=\Psi(x,y),(x,y)\in \Omega \space\space\space\space(1),\\ u(0,y,t)=g_{1}(y,t),u(a,y,t)=g_{2}(y,t),0\leqslant y \leqslant b,0<t\leqslant T,\\ u(x,0,t)=g_{3}(x,t),u(x,b,t)=g_{4}(x,t),0\leqslant x \leqslant a,0<t \leqslant T \end{matrix}\right.

        为了保证连续性,公式(1)中相关函数满足:

g_{1}(0,t)=g_{3}(0,t),g_{1}(b,t)=g_{4}(b,t)\\g_{2}(0,t)=g_{3}(a,t),g_{2}(b,t)=g_{4}(a,t)

\varphi(0,y)=g_{1}(y,0),\varphi(a,y)=g_{2}(y,0)\\\varphi(x,0)=g_{3}(x,0),\varphi(x,b)=g_{4}(x,0)

二、理论推导

2.1 显差分格式

        第一步:网格剖分。在三维长方体空间(二维平面加一维时间轴)进行网格剖分,将区域[0,a]等分m份,区域[0,b]等分n份,区域[0,T]等分l份,即

x_{i}=i\cdot\Delta x=\frac{ia}{m},0\leqslant i\leqslant m,y_{j}=j\cdot\Delta y=\frac{jb}{n},0\leqslant j\leqslant n,t_{k}=k\cdot\Delta t=\frac{kT}{l},0\leqslant k\leqslant l

得到网格节点坐标为(x_{i},y_{j},t_{k})。我们利用数值方法获得精确解u(x,y,t)在网格节点(x_{i},y_{j},t_{k})处的近似值,即数值解u^{k}_{i,j}

        第二步:弱化原方程,仅在离散点处成立。即

\left\{\begin{matrix} \frac{\partial^{2}u}{\partial t^{2}}|_{(x_{i},y_{j},t_{k})}-(\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}}{\partial y^{2}})|_{(x_{i},y_{j},t_{k})}=f(x_{i},y_{j},t_{k}),0\leqslant i\leqslant m,0\leqslant j\leqslant n,0<k\leqslant l,\\ u(x_{i},y_{j},t_{0})=\varphi(x_{i},y_{j}),\frac{\partial u}{\partial t}(x_{i},y_{j},t_{0})=\Psi(x_{i},y_{j}),0\leqslant i\leqslant m,0\leqslant j\leqslant n,\space\space(2)\\ u(x_{0},y_{j},t_{k})=g_{1}(y_{j},t_{k}),u(x_{m},y_{j},t_{k})=g_{2}(y_{j},t_{k}),0\leqslant j\leqslant n,0<k\leqslant l,\\ u(x_{i},y_{0},t_{k})=g_{3}(x_{i},t_{k}),u(x_{i},y_{n},t_{k})=g_{4}(x_{i},t_{k}),0\leqslant i\leqslant m,0<k\leqslant l \end{matrix}\right.

        第三步:差商代替微商处理偏导数。方程组二阶偏导数使用中心差商来近似,即

\frac{\partial^{2}u}{\partial t^{2}}|_{(x_{i},y_{j},t_{k})}\approx\frac{u(x_{i},y_{j},t_{k-1})-2u(x_{i},y_{j},t_{k})+u(x_{i},y_{j},t_{k+1})}{\Delta t^{2}}

\frac{\partial^{2}u}{\partial x^{2}}|_{(x_{i},y_{j},t_{k})}\approx\frac{u(x_{i-1},y_{j},t_{k})-2u(x_{i},y_{j},t_{k})+u(x_{i+1},y_{j},t_{k})}{\Delta x^{2}}

\frac{\partial^{2}u}{\partial y^{2}}|_{(x_{i},y_{j},t_{k})}\approx\frac{u(x_{i},y_{j-1},t_{k})-2u(x_{i},y_{j},t_{k})+u(x_{i},y_{j+1},t_{k})}{\Delta y^{2}}

对于初始条件中的一阶偏导数,也采用二阶中心差商,即

\frac{\partial u}{\partial t}|_{(x_{i},y_{j},t_{k})}\approx\frac{u(x_{i},y_{j},t_{1})-u(x_{i},y_{j},t_{-1})}{2\Delta t}

        将差商公式带入公式(2),用数值解代替精确解并忽略高阶项,可得对应的数值格式为

\left\{\begin{matrix} \frac{u^{k-1}_{i,j}-2u^{k}_{i,j}+u^{k+1}_{i,j}}{\Delta t^{2}}-(\frac{u^{k}_{i-1,j}-2u^{k}_{i,j}+u^{k}_{i+1,j}}{\Delta x^{2}}+\frac{u^{k}_{i,j-1}-2u^{k}_{i,j}+u^{k}_{i,j+1}}{\Delta y^{2}})=f(x_{i},y_{j},t_{k}),1\leqslant i\leqslant m-1,1\leqslant j\leqslant n-1,1\leqslant k\leqslant l-1,\\ u^{0}_{i,j}=\varphi(x_{i},y_{j}),u^{1}_{i,j}=u^{-1}_{i,j}+2\Delta t\Psi(x_{i},y_{j}),0\leqslant i\leqslant m,0\leqslant j\leqslant n,\space\space(3)\\ u^{k}_{0,j}=g_{1}(y_{j},t_{k}),u^{k}_{m,j}=g_{2}(y_{j},t_{k}),0\leqslant j\leqslant n,0<k\leqslant l,\\ u^{k}_{i,0}=g_{3}(x_{i},t_{k}),u^{k}_{i,n}=g_{4}(x_{i},t_{k}),0\leqslant i\leqslant m,0<k\leqslant l \end{matrix}\right.

        为消去虚拟边界项u^{-1}_{i,j},让公式(3)中第1式对k=0成立,可得

\frac{u^{-1}_{i,j}-2u^{0}_{i,j}+u^{1}_{i,j}}{\Delta t^{2}}-(\frac{u^{0}_{i-1,j}-2u^{0}_{i,j}+u^{0}_{i+1,j}}{\Delta x^{2}}+\frac{u^{0}_{i,j-1}-2u^{0}_{i,j}+u^{0}_{i,j+1}}{\Delta y^{2}})=f(x_{i},y_{j},t_{0})

        再结合公式(3)中第3式中越界的u^{-1}_{i,j}=u^{1}_{i,j}-2\Delta t\Psi(x_{i},y_{j}),消去虚拟边界项u^{-1}_{i,j},可得

2u^{1}_{i,j}-2\Delta t\Psi(x_{i},y_{j})-2u^{0}_{i,j}=\Delta t^{2}(\frac{u^{0}_{i-1,j}-2u^{0}_{i,j}+u^{0}_{i+1,j}}{\Delta x^{2}}+\frac{u^{0}_{i,j-1}-2u^{0}_{i,j}+u^{0}_{i,j+1}}{\Delta y^{2}}+f(x_{i},y_{j},t_{0}))

r_{1}=\Delta t^{2}/\Delta x^{2}, r_{2}=\Delta t^{2}/\Delta y^{2},整理上式可得

u^{1}_{i,j}=\Delta t\Psi(x_{i},y_{j})+(r_{1}(u^{0}_{i-1,j}+u^{0}_{i+1,j})+r_{2}(u^{0}_{i,j-1}+u^{0}_{i,j+1})+f(x_{i},y_{j},t_{0})\Delta t^{2})/2+(1-r_{1}-r_{2})u^{0}_{i,j} \space\space\space\space(4)

        由公式(4)知,公式(3)即为三层显格式:

\left\{\begin{matrix} u^{k+1}_{i,j}=r_{1}(u^{k}_{i-1,j}+u^{k}_{i+1,j})+r_{2}(u^{0}_{i,j-1}+u^{k}_{i,j+1})+2(1-r_{1}-r_{2})u^{k}_{i,j}-u^{k-1}_{i,j}=f(x_{i},y_{j},t_{k})\Delta t^{2},1\leqslant i\leqslant m-1,1\leqslant j\leqslant n-1,1\leqslant k\leqslant l-1,\\ u^{0}_{i,j}=\varphi(x_{i},y_{j}),0\leqslant i\leqslant m,0\leqslant j\leqslant n,\\ u^{1}_{i,j}=\Delta t\Psi(x_{i},y_{j})+(r_{1}(u^{0}_{i-1,j}+u^{0}_{i+1,j})+r_{2}(u^{0}_{i,j-1}+u^{0}_{i,j+1})+f(x_{i},y_{j},t_{0})\Delta t^{2})/2+(1-r_{1}-r_{2})u^{0}_{i,j},1\leqslant i\leqslant m-1,1\leqslant j\leqslant n-1,\\ u^{k}_{0,j}=g_{1}(y_{j},t_{k}),u^{k}_{m,j}=g_{2}(y_{j},t_{k}),0\leqslant j\leqslant n,0<k\leqslant l,\\ u^{k}_{i,0}=g_{3}(x_{i},t_{k}),u^{k}_{i,n}=g_{4}(x_{i},t_{k}),0\leqslant i\leqslant m,0<k\leqslant l \end{matrix}\right.

该格式的局部截断误差为O(\Delta t^{2}+\Delta x^{2}+\Delta y^{2}),考虑到其稳定性限制,将其改为隐格式。

2.2 交替方向隐格式

        由于

\frac{\partial ^{2}u}{\partial x^{2}}|_{(x_{i,}y_{j},t_{k})}\approx\frac{1}{2}(\frac{\partial ^{2}u}{\partial x^{2}}|_{(x_{i},y_{j},t_{k-1})}+\frac{\partial^{2}u}{\partial x^{2}}|_{(x_{i},y_{j},t_{k+1})})

\frac{\partial ^{2}u}{\partial y^{2}}|_{(x_{i,}y_{j},t_{k})}\approx\frac{1}{2}(\frac{\partial ^{2}u}{\partial y^{2}}|_{(x_{i},y_{j},t_{k-1})}+\frac{\partial^{2}u}{\partial y^{2}}|_{(x_{i},y_{j},t_{k+1})})

由上面两式可将原来在点(x_{i},y_{j},t_{k})处的连续方程

\frac{\partial^{2}u}{\partial t^{2}}|_{(x_{i},y_{j},t_{k})}-(\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}})|_{(x_{i},y_{j},t_{k})}=f(x_{i},y_{j},t_{k})

改写为

\frac{\partial^{2}u}{\partial t^{2}}|_{(x_{i},y_{j},t_{k})}-\frac{1}{2}((\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}})|_{(x_{i},y_{j},t_{k-1})}+(\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}})|_{(x_{i},y_{j},t_{k+1})})=f(x_{i},y_{j},t_{k})+O(\Delta t^{2})

        再用中心差商处理二阶偏导数,将数值解代替精确解并忽略高阶项,可得数值隐格式为:

\frac{u^{k-1}_{i,j}-2u^{k}_{i,j}+u^{k+1}_{i,j}}{\Delta t^{2}}-\frac{1}{2}(\frac{u^{k-1}_{i-1,j}-2u^{k-1}_{i,j}+u^{k-1}_{i+1,j}}{\Delta x^{2}}+\frac{u^{k+1}_{i-1,j}-2u^{k+1}_{i,j}+u^{k+1}_{i+1,j}}{\Delta x^{2}})-\frac{1}{2}(\frac{u^{k-1}_{i,j-1}-2u^{k-1}_{i,j}+u^{k-1}_{i,j+1}}{\Delta y^{2}}+\frac{u^{k+1}_{i,j-1}-2u^{k+1}_{i,j}+u^{k+1}_{i,j+1}}{\Delta y^{2}})=f(x_{i},y_{j},t_{k}) \space\space(5)

        公式(5)在求解上存在困难,故参照二维抛物型偏微分方程初边值问题的交替方向隐格式设计思想,先将公式(5)写成算子的形式,可改写为:

\frac{u^{k-1}_{i,j}-2u^{k}_{i,j}+u^{k+1}_{i,j}}{\Delta t^{2}}-\frac{1}{2}(\frac{\delta^{2}_{x}u^{k-1}_{i,j}}{\Delta x^{2}}+\frac{\delta^{2}_{x}u^{k+1}_{i,j}}{\Delta x^{2}})-\frac{1}{2}(\frac{\delta^{2}_{y}u^{k-1}_{i,j}}{\Delta y^{2}}+\frac{\delta^{2}_{y}u^{k+1}_{i,j}}{\Delta y^{2}})=f(x_{i},y_{j},t_{k})

也就是          u^{k+1}_{i,j}-2^{k}_{i,j}+u^{k+1}_{i,j}=\frac{r_{1}}{2}(\delta^{2}_{x}u^{k-1}_{i,j}+\delta^{2}_{x}u^{k+1}_{i,j})+\frac{r_{2}}{2}(\delta^{2}_{y}u^{k-1}_{i,j}+\delta^{2}_{y}u^{k+1}_{i,j})+f(x_{i},y_{j},t_{k})\Delta t^{2}

(1-\frac{r_{1}}{2}\delta^{2}_{x}-\frac{r_{2}}{2}\delta^{2}_{y})u^{k+1}_{i,j}=(\frac{r_{1}}{2}\delta^{2}_{x}+\frac{r_{2}}{2}\delta^{2}_{y}-1)u^{k-1}_{i,j}+2u^{k}_{i,j}+f(x_{i},y_{j},t_{k})\Delta t^{2}

        对上式左端添加一个辅助项,即

(1-\frac{r_{1}}{2}\delta^{2}_{x}-\frac{r_{2}}{2}\delta^{2}_{y}+\frac{r_{1}r_{2}}{4}\delta^{2}_{x}\delta^{2}_{y})u^{k+1}_{i,j}=(\frac{r_{1}}{2}\delta^{2}_{x}+\frac{r_{2}}{2}\delta^{2}_{y}-1)u^{k-1}_{i,j}+2u^{k}_{i,j}+f(x_{i},y_{j},t_{k})\Delta t^{2} \space\space(6)

        对公式(6)进行分解,可以写为

(1-\frac{r_{1}}{2}\delta^{2}_{x})(1-\frac{r_{2}}{2}\delta^{2}_{y})u^{k+1}_{i,j}=(\frac{r_{1}}{2}\delta^{2}_{x}+\frac{r_{2}}{2}\delta^{2}_{y}-1)u^{k-1}_{i,j}+2u^{k}_{i,j}+f(x_{i},y_{j},t_{k})\Delta t^{2} \space\space(6)

        通过引入中间变量V_{i,j},得

\left\{\begin{matrix} (1-\frac{r_{1}}{2}\delta^{2}_{x})V_{i,j}=(\frac{r_{1}}{2}\delta^{2}_{x}+\frac{r_{2}}{2}\delta^{2}_{y}-1)u^{k-1}_{i,j}+2u^{k}_{i,j}+f(x_{i},y_{j},t_{k})\\ (1-\frac{r_{2}}{2}\delta^{2}_{y})u^{k+1}_{i,j}=V_{i,j} \space\space(7)\end{matrix}\right.

        上式是一个x方向的三对角线性方程组,计算第1式中的V_{i,j}时需要用到V_{0,j}V_{m,j},这可以直接从第2式中解得,即

V_{0,j}=(1-\frac{r_{2}}{2}\delta^{2}_{y})u^{k+1}_{0,j},V_{m,j}=(1-\frac{r_{2}}{2}\delta^{2}_{y})u^{k+1}_{m,j} \space\space(8)

        我们可以直接将公式(7)、(8)写成差分形式,即

\left\{\begin{matrix} -\frac{r_{1}}{2}V_{i-1,j}+(1+r_{1})V_{i,j}-\frac{r_{1}}{2}V_{i+1,j}=2u^{k}_{i,j}+f(x_{i},y_{j},t_{k})\Delta t^{2}+\\\frac{r_{1}}{2}(u^{k-1}_{i-1,j}+u^{k-1}_{i+1,j})+\frac{r_{2}}{2}(u^{k-1}_{i,j-1}+u^{k-1}_{i,j+1})-(1+r_{1}+r_{2})u^{k-1}_{i,j}\\ -\frac{r_{2}}{2}u^{k+1}_{i,j-1}+(1+r_{2})u^{k+1}_{i,j}-\frac{r_{2}}{2}u^{k+1}_{i,j+1}=V_{i,j} \space\space(9)\end{matrix}\right.

其中,                            V_{0,j}=-\frac{r_{2}}{2}(u^{k+1}_{0,j-1}+u^{k+1}_{0,j+1})+(1+r_{2})u^{k+1}_{0,j}

V_{m,j}=-\frac{r_{2}}{2}(u^{k+1}_{m,j-1}+u^{k+1}_{m,j+1})+(1+r_{2})u^{k+1}_{m,j}

        假设第k和k-1个时间层的信息已知,这样在公式(9)的第1式中,先固定某个j,这时公式(9)的第1式就是一个m-1阶线性方程组,其系数矩阵是三对角矩阵,可以用追赶法求解。每解一个方程组就可以得到一个在[0,a]x[0,b]矩形区域内行网格点(即行节点(x_{i},y_{j}))上的中间量V的信息,所以要得到矩形区域内所有的网格点信息,需要求解n-1个这样的线性系统。同样,公式(9)的第2式是一个y方向的三对角线性方程组时,求解时要先固定某个i,这时公式(9)的第2式是一个n-1阶线性方程组,系数矩阵是三对角矩阵,可以用追赶法求解,每解一个方程组得到第k+1个时间层上[0,a]x[0,b]矩形区域内列网格点(即(x_{i},y_{j},t_{k+1}))上的数值解的信息,所以要得到第k+1时间层上所有网格点的信息,需要求解m-1个这样的线性系统。这样就实现了x方向和y方向交替求解的隐格式。

三、算例实现

        用交替方向隐格式求解下列二维双曲型方程初边值问题:

\left\{\begin{matrix} \frac{\partial^{2}u}{\partial t^{2}}-(\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial y^{2}})=\frac{4t}{(1+x^{2}+y^{2})^{2}}(1-\frac{2(x^{2}+y^{2})}{1+x^{2}+y^{2}}),0<x,y<1,0<t\leqslant 1,\\ u(x,y,0)=0,\frac{\partial u}{\partial t}(x,y,0)=\frac{1}{1+x^{2}+y^{2}},0\leqslant x,y\leqslant 1,\\ u(0,y,t)=\frac{t}{1+y^{2}},u(1,y,t)=\frac{t}{2+y^{2}},0\leqslant y\leqslant 1,0<t\leqslant 1,\\ u(x,0,t)=\frac{t}{1+x^{2}},u(x,1,t)=\frac{t}{2+x^{2}},0<x<1,0<t\leqslant 1. \end{matrix}\right.

已知精确解为u(x,y,t)=\frac{t}{1+x^{2}+y^{2}}。分别取步长为\Delta x=1/10,\Delta y=1/20,\Delta t=1/40\Delta x=1/20,\Delta y=1/40,\Delta t=1/80,给出在节点(0.2i,0.25j,1.00),i=1,2,3,4,j=1,2,3处的数值解和误差。

代码如下:

#include<cmath>
#include<stdlib.h>
#include<stdio.h>int main(int argc, char* argv[])
{int i, j, k, m, n, L, gap_i, gap_j;double a, b, T, r1, r2, dx, dy, dt, *x, *y, *t, ***u, **v;double *a1, *b1, *c1, *d1, *a2, *b2, *c2, *d2, *ans, temp;double f(double x, double y, double t);double phi(double x, double y);double psi(double x, double y);double g1(double y, double t);double g2(double y, double t);double g3(double x, double t);double g4(double x, double t);double *chase_algorithm(double *a, double *b, double *c, double *d, int n);double exact(double x, double y, double t);a = 1.0;b = 1.0;T = 1.0;m = 10;n = 20;L = 40;dx = a / m;dy = b / n;dt = T / L;temp = dt / dx;r1 = temp*temp;temp = dt / dy;r2 = temp*temp;printf("m=%d, n=%d, L=%d.\n", m, n, L);printf("r1=%.4f,   r2=%.4f.\n", r1, r2);x = (double *)malloc(sizeof(double)*(m + 1));for (i = 0; i <= m; i++)x[i] = i*dx;y = (double *)malloc(sizeof(double)*(n + 1));for (j = 0; j <= n; j++)y[j] = j*dy;t = (double *)malloc(sizeof(double)*(L + 1));for (k = 0; k <= L; k++)t[k] = k*dt;u = (double ***)malloc(sizeof(double **)*((m + 1)*(n + 1)*(L + 1)));v = (double **)malloc(sizeof(double *)*((m + 1)*(n + 1)));for (i = 0; i <= m; i++){u[i] = (double **)malloc(sizeof(double *)*((n + 1)*(L + 1)));v[i] = (double *)malloc(sizeof(double)*(n + 1));}for (i = 0; i <= m; i++){for (j = 0; j <= n; j++)u[i][j] = (double *)malloc(sizeof(double)*(L + 1));}for (i = 0; i <= m; i++){for (j = 0; j <= n; j++){u[i][j][0] = phi(x[i], y[j]);}}for (i = 1; i< m; i++){for (j = 1; j < n; j++){u[i][j][1] = u[i][j][0] + dt*psi(x[i], y[j]) + (r1*(u[i - 1][j][0] - 2 * u[i][j][0] + u[i + 1][j][0]) + r2*(u[i][j - 1][0] - 2 * u[i][j][0] + u[i][j + 1][0]) + dt*dt*f(x[i], y[j], t[0])) / 2.0;}}for (k = 1; k <= L; k++){for (j = 0; j <= n; j++){u[0][j][k] = g1(y[j], t[k]);u[m][j][k] = g2(y[j], t[k]);}for (i = 1; i <= m - 1; i++){u[i][0][k] = g3(x[i], t[k]);u[i][n][k] = g4(x[i], t[k]);}}a1 = (double *)malloc(sizeof(double)*(m - 1));b1 = (double *)malloc(sizeof(double)*(m - 1));c1 = (double *)malloc(sizeof(double)*(m - 1));d1 = (double *)malloc(sizeof(double)*(m - 1));for (i = 0; i< m - 1; i++){a1[i] = -r1 / 2.0;b1[i] = 1.0 + r1;c1[i] = a1[i];}a2 = (double *)malloc(sizeof(double)*(n - 1));b2 = (double *)malloc(sizeof(double)*(n - 1));c2 = (double *)malloc(sizeof(double)*(n - 1));d2 = (double *)malloc(sizeof(double)*(n - 1));for (j = 0; j < n - 1; j++){a2[j] = -r2 / 2.0;b2[j] = 1.0 + r2;c2[j] = a2[j];}for (k = 1; k < L; k++){for (j = 1; j <= n - 1; j++){for (i = 1; i <= m - 1; i++){d1[i-1]=2 * u[i][j][k] + dt*dt*f(x[i], y[j], t[k])+r1*(u[i-1][j][k-1]+u[i+1][j][k-1])/2.0+r2*(u[i][j-1][k-1]+u[i][j+1][k-1])/2.0-(1+r1+r2)*u[i][j][k-1];}v[0][j] = -r2*(u[0][j - 1][k + 1] + u[0][j + 1][k + 1]) / 2.0 + (1 + r2)*u[0][j][k + 1];d1[0] = d1[0] + r1*v[0][j] / 2.0;v[m][j] = -r2*(u[m][j - 1][k + 1] + u[m][j + 1][k + 1]) / 2.0 + (1 + r2)*u[m][j][k + 1];d1[m - 2] = d1[m - 2] + r1*v[m][j] / 2.0;ans = chase_algorithm(a1, b1, c1, d1, m - 1);for (i = 1; i <= m - 1; i++)v[i][j] = ans[i - 1];free(ans);}for (i = 1; i <= m - 1; i++){for (j = 1; j <= n - 1; j++)d2[j - 1] = v[i][j];d2[0] = d2[0] + r2*u[i][0][k + 1] / 2.0;d2[n - 2] = d2[n - 2] + r2*u[i][n][k + 1] / 2.0;ans = chase_algorithm(a2, b2, c2, d2, n - 1);for (j = 1; j <= n - 1; j++)u[i][j][k + 1] = ans[j - 1];free(ans);}}k = L ;gap_i = m / 5;gap_j = n / 4;for (i = gap_i; i <= m - 1; i = i + gap_i){for (j = gap_j; j <= n - 1; j = j + gap_j){temp = fabs(exact(x[i], y[j], t[k]) - u[i][j][k]);printf("(%.2f, %.2f, %.2f)   y=%f,  err=%.4e\n", x[i], y[j], t[k], u[i][j][k], temp);}}free(x); free(y); free(t);free(a1); free(b1); free(c1); free(d1);free(a2); free(b2); free(c2); free(d2);for (i = 0; i <= m; i++){for (j = 0; j <= n; j++)free(u[i][j]);}for (i = 0; i <= m; i++){free(u[i]);}free(u);return 0;
}double f(double x, double y, double t)
{double z,temp;z = x*x + y*y;temp = 1 - 2 * z / (1 + z);return 4*t*temp/((1+z)*(1+z));
}
double phi(double x, double y)
{return 0.0;
}
double psi(double x, double y)
{return 1.0/(1+x*x+y*y);
}
double g1(double y, double t)
{return t/(1+y*y);
}
double g2(double y, double t)
{return t/(2+y*y);
}
double g3(double x, double t)
{return t/(1+x*x);
}
double g4(double x, double t)
{return t/(2+x*x);
}
double *chase_algorithm(double *a, double *b, double *c, double *d, int n)
{int i;double * ans, *g, *w, p;ans = (double *)malloc(sizeof(double)*n);g = (double *)malloc(sizeof(double)*n);w = (double *)malloc(sizeof(double)*n);g[0] = d[0] / b[0];w[0] = c[0] / b[0];for (i = 1; i<n; i++){p = b[i] - a[i] * w[i - 1];g[i] = (d[i] - a[i] * g[i - 1]) / p;w[i] = c[i] / p;}ans[n - 1] = g[n - 1];i = n - 2;do{ans[i] = g[i] - w[i] * ans[i + 1];i = i - 1;} while (i >= 0);free(g);free(w);return ans;
}
double exact(double x, double y, double t)
{return t/(1+x*x+y*y);
}

\Delta x=1/10,\Delta y=1/20,\Delta t=1/40时,计算结果如下:

m=10, n=20, L=40.
r1=0.0625,   r2=0.2500.
(0.20, 0.25, 1.00)   y=0.907133,  err=1.0342e-04
(0.20, 0.50, 1.00)   y=0.775225,  err=3.0722e-05
(0.20, 0.75, 1.00)   y=0.624030,  err=4.5462e-06
(0.40, 0.25, 1.00)   y=0.817921,  err=7.4625e-05
(0.40, 0.50, 1.00)   y=0.709103,  err=1.1636e-04
(0.40, 0.75, 1.00)   y=0.580475,  err=7.6977e-05
(0.60, 0.25, 1.00)   y=0.702808,  err=1.7953e-04
(0.60, 0.50, 1.00)   y=0.620917,  err=2.0139e-04
(0.60, 0.75, 1.00)   y=0.520023,  err=1.3270e-04
(0.80, 0.25, 1.00)   y=0.587236,  err=1.3524e-04
(0.80, 0.50, 1.00)   y=0.528946,  err=1.5470e-04
(0.80, 0.75, 1.00)   y=0.453919,  err=1.1003e-04

\Delta x=1/20,\Delta y=1/40,\Delta t=1/80时,计算结果如下:

m=20, n=40, L=80.
r1=0.0625,   r2=0.2500.
(0.20, 0.25, 1.00)   y=0.907056,  err=2.6327e-05
(0.20, 0.50, 1.00)   y=0.775202,  err=7.7268e-06
(0.20, 0.75, 1.00)   y=0.624026,  err=9.2774e-07
(0.40, 0.25, 1.00)   y=0.817978,  err=1.8204e-05
(0.40, 0.50, 1.00)   y=0.709191,  err=2.8662e-05
(0.40, 0.75, 1.00)   y=0.580532,  err=1.9294e-05
(0.60, 0.25, 1.00)   y=0.702943,  err=4.4834e-05
(0.60, 0.50, 1.00)   y=0.621068,  err=5.0469e-05
(0.60, 0.75, 1.00)   y=0.520123,  err=3.2859e-05
(0.80, 0.25, 1.00)   y=0.587338,  err=3.3376e-05
(0.80, 0.50, 1.00)   y=0.529063,  err=3.7709e-05
(0.80, 0.75, 1.00)   y=0.454003,  err=2.6860e-05

四、结论

        从计算结果可以看出,将时间、空间步长同时减小1/2,误差减小1/4,可见该数值格式是二阶的。

这篇关于偏微分方程算法之二阶双曲型方程交替方向隐格式的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/933497

相关文章

不懂推荐算法也能设计推荐系统

本文以商业化应用推荐为例,告诉我们不懂推荐算法的产品,也能从产品侧出发, 设计出一款不错的推荐系统。 相信很多新手产品,看到算法二字,多是懵圈的。 什么排序算法、最短路径等都是相对传统的算法(注:传统是指科班出身的产品都会接触过)。但对于推荐算法,多数产品对着网上搜到的资源,都会无从下手。特别当某些推荐算法 和 “AI”扯上关系后,更是加大了理解的难度。 但,不了解推荐算法,就无法做推荐系

康拓展开(hash算法中会用到)

康拓展开是一个全排列到一个自然数的双射(也就是某个全排列与某个自然数一一对应) 公式: X=a[n]*(n-1)!+a[n-1]*(n-2)!+...+a[i]*(i-1)!+...+a[1]*0! 其中,a[i]为整数,并且0<=a[i]<i,1<=i<=n。(a[i]在不同应用中的含义不同); 典型应用: 计算当前排列在所有由小到大全排列中的顺序,也就是说求当前排列是第

csu 1446 Problem J Modified LCS (扩展欧几里得算法的简单应用)

这是一道扩展欧几里得算法的简单应用题,这题是在湖南多校训练赛中队友ac的一道题,在比赛之后请教了队友,然后自己把它a掉 这也是自己独自做扩展欧几里得算法的题目 题意:把题意转变下就变成了:求d1*x - d2*y = f2 - f1的解,很明显用exgcd来解 下面介绍一下exgcd的一些知识点:求ax + by = c的解 一、首先求ax + by = gcd(a,b)的解 这个

综合安防管理平台LntonAIServer视频监控汇聚抖动检测算法优势

LntonAIServer视频质量诊断功能中的抖动检测是一个专门针对视频稳定性进行分析的功能。抖动通常是指视频帧之间的不必要运动,这种运动可能是由于摄像机的移动、传输中的错误或编解码问题导致的。抖动检测对于确保视频内容的平滑性和观看体验至关重要。 优势 1. 提高图像质量 - 清晰度提升:减少抖动,提高图像的清晰度和细节表现力,使得监控画面更加真实可信。 - 细节增强:在低光条件下,抖

【数据结构】——原来排序算法搞懂这些就行,轻松拿捏

前言:快速排序的实现最重要的是找基准值,下面让我们来了解如何实现找基准值 基准值的注释:在快排的过程中,每一次我们要取一个元素作为枢纽值,以这个数字来将序列划分为两部分。 在此我们采用三数取中法,也就是取左端、中间、右端三个数,然后进行排序,将中间数作为枢纽值。 快速排序实现主框架: //快速排序 void QuickSort(int* arr, int left, int rig

poj 3974 and hdu 3068 最长回文串的O(n)解法(Manacher算法)

求一段字符串中的最长回文串。 因为数据量比较大,用原来的O(n^2)会爆。 小白上的O(n^2)解法代码:TLE啦~ #include<stdio.h>#include<string.h>const int Maxn = 1000000;char s[Maxn];int main(){char e[] = {"END"};while(scanf("%s", s) != EO

秋招最新大模型算法面试,熬夜都要肝完它

💥大家在面试大模型LLM这个板块的时候,不知道面试完会不会复盘、总结,做笔记的习惯,这份大模型算法岗面试八股笔记也帮助不少人拿到过offer ✨对于面试大模型算法工程师会有一定的帮助,都附有完整答案,熬夜也要看完,祝大家一臂之力 这份《大模型算法工程师面试题》已经上传CSDN,还有完整版的大模型 AI 学习资料,朋友们如果需要可以微信扫描下方CSDN官方认证二维码免费领取【保证100%免费

【生成模型系列(初级)】嵌入(Embedding)方程——自然语言处理的数学灵魂【通俗理解】

【通俗理解】嵌入(Embedding)方程——自然语言处理的数学灵魂 关键词提炼 #嵌入方程 #自然语言处理 #词向量 #机器学习 #神经网络 #向量空间模型 #Siri #Google翻译 #AlexNet 第一节:嵌入方程的类比与核心概念【尽可能通俗】 嵌入方程可以被看作是自然语言处理中的“翻译机”,它将文本中的单词或短语转换成计算机能够理解的数学形式,即向量。 正如翻译机将一种语言

dp算法练习题【8】

不同二叉搜索树 96. 不同的二叉搜索树 给你一个整数 n ,求恰由 n 个节点组成且节点值从 1 到 n 互不相同的 二叉搜索树 有多少种?返回满足题意的二叉搜索树的种数。 示例 1: 输入:n = 3输出:5 示例 2: 输入:n = 1输出:1 class Solution {public int numTrees(int n) {int[] dp = new int

Codeforces Round #240 (Div. 2) E分治算法探究1

Codeforces Round #240 (Div. 2) E  http://codeforces.com/contest/415/problem/E 2^n个数,每次操作将其分成2^q份,对于每一份内部的数进行翻转(逆序),每次操作完后输出操作后新序列的逆序对数。 图一:  划分子问题。 图二: 分而治之,=>  合并 。 图三: 回溯: