偏微分方程算法之九点紧差分法

2024-04-30 22:12

本文主要是介绍偏微分方程算法之九点紧差分法,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

目录

一、研究目标

二、理论推导

三、算例实现

四、结论


一、研究目标

        我们已经在专栏中介绍了椭圆型偏微分方程的五点菱形差分格式,这里我们继续以该方法为背景,探讨如何提高五点法的精度,即从二阶精度提升到四阶精度。

        研究目标现继续以矩形区域\Omega=[(x,y)|a\leqslant x\leqslant b,c\leqslant y\leqslant d]内的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)\\ u(x,y)=\varphi(x,y),(x,y)\in\partial \Omega=\Gamma \end{matrix}\right.

二、理论推导

        提高差分格式精度的一种有效方法是紧差分,即引入中间函数,对二阶偏导数进行更精细的逼近。

        与五点菱形差分法类似,首先进行矩形区域的剖分,然后将原方程弱化,使其仅在离散节点处成立。在利用差商代替微商之前,引入中间函数,令

\frac{\partial ^{2}u}{\partial x^{2}}=v(x,y),\frac{\partial^{2}u}{\partial y^{2}}=w(x,y)

原方程在离散节点处满足

-(v+w)|_{(x_{i},y_{j})}=f(x_{i},y_{j})\; (2)

又有,当 max_{(x,y)\in\Omega}[|\frac{\partial^{6}u(x,y)}{\partial x^{6}}|,|\frac{\partial^{6}u(x,y)}{\partial y^{6}}|]有界时,

\frac{\partial^{2}u}{\partial x^{2}}|_{(x_{i},y_{j})}=\frac{u(x_{i-1},y_{j})-2u(x_{i},y_{j})+u(x_{i+1},y_{j})}{\Delta x^{2}}-\frac{\Delta x^{2}}{12}\frac{\partial^{4}u}{\partial x^{4}}|_{(x_{i},y_{j})}+C_{1}(\Delta x)^{4}

从而有

v(x_{i},y_{j})=\frac{\partial^{2}u}{\partial x^{2}}|_{(x_{i},y_{j})}=\frac{\delta ^{2}_{x}u(x_{i},y_{j})}{\Delta x^{2}}-\frac{\Delta x^{2}}{12}\frac{\partial^{2}v}{\partial x^{2}}|_{(x_{i},y_{j})}+C_{1}(\Delta x)^{4}

=\frac{\delta^{2}_{x}u(x_{i},y_{j})}{\Delta x^{2}}-\frac{\Delta x^{2}}{12}(\frac{v(x_{i-1},y_{j})-2v(x_{i},y_{j})+v(x_{i+1},y_{j})}{\Delta x^{2}})+C_{2}(\Delta x)^{4}

也就是

v(x_{i},y_{j})=\frac{\delta^{2}_{x}u(x_{i},y_{j})}{\Delta x^{2}}-\frac{1}{12}(v(x_{i-1},y_{j})+10v(x_{i},y_{j})+v(x_{i+1},y_{j}))+C_{2}(\Delta x)^{4}\; (3)

由二维抛物型偏微分方程紧交替方向隐格式方法介绍中算子公式

\varepsilon ^{2}_{x}v(x_{i},y,t)=\frac{1}{12}(v(x_{i-1},y,t)+10v(x_{i},y,t)+v(x_{i+1},y,t))

可知公式(3)为

\epsilon^{2}_{x}v(x_{i},y_{j})=\frac{\delta^{2}_{x}u(x_{i},y_{j})}{\Delta x^{2}}+C_{2}(\Delta x)^{4}\; (4)

同理可得

\epsilon^{2}_{y}w(x_{i},y_{j})=\frac{\delta^{2}_{y}u(x_{i},y_{j})}{\Delta y^{2}}+C_{3}(\Delta y)^{4}\; (5)

        将算子\epsilon^{2}_{x}\epsilon^{2}_{y}作用到公式(2)上,有

-\epsilon^{2}_{x}\epsilon^{2}_{y}(v+w)|_{(x_{i},y_{j})}=\epsilon^{2}_{x}\epsilon^{2}_{y}f(x_{i},y_{j})\; (6)

        然后将公式(4)、(5)一起代入公式(6),可得

-\frac{\epsilon^{2}_{y}\delta^{2}_{x}u(x_{i},y_{j})}{\Delta x^{2}}-\frac{\epsilon^{2}_{x}\delta^{2}_{y}u(x_{i},y_{j})}{\Delta y^{2}}=\epsilon^{2}_{x}\epsilon^{2}_{y}f(x_{i},y_{j})+C_{2}(\Delta x)^{4}+C_{3}(\Delta y)^{4}

        用数值解代替精确解并忽略高阶项,可得紧差分格式:

\left\{\begin{matrix} -\frac{\epsilon^{2}_{y}\delta^{2}_{x}u_{i,j}}{\Delta x^{2}}-\frac{\epsilon^{2}_{x}\delta^{2}_{y}u_{i,j}}{\Delta y^{2}}=\epsilon^{2}_{x}\epsilon^{2}_{y}f(x_{i},y_{j}),1\leqslant i\leqslant m-1,1\leqslant j\leqslant n-1,\; (7)\\ u_{s,t}=\varphi(x_{s},y_{t}),s=0,m,0\leqslant t\leqslant n;t=0,n,0\leqslant s\leqslant m \end{matrix}\right.

该格式的局部截断误差为O(\Delta x^{4}+\Delta y^{4}),将公式(7)中第1式整理化简,利用记号

\frac{1}{\Delta x^{2}}=\beta,\frac{1}{\Delta y^{2}}=\gamma

可得

(\beta+\gamma)u_{i-1,j-1}+(10\gamma-2\beta)u_{i,j-1}+(\beta+\gamma)u_{i+1,j-1}+(10\beta-2\gamma)u_{i-1,j}-20(\beta+\gamma)u_{i,j}+(10\beta-2\gamma)u_{i+1,j}+(\beta+\gamma )u_{i-1,j+1}+(10\gamma-2\beta)u_{i,j+1}+(\beta+\gamma)u_{i+1,j+1}=-12\epsilon^{2}_{x}\epsilon^{2}_{y}f(x_{i},y_{j})

        该紧格式每一步计算都涉及9个点,所以公式(7)称为九点紧差分格式。

        该格式无法写成线性方程组\mathbf{A}x=\mathbf{b}的形式,只能写成

\begin{pmatrix} \eta_{2} & \xi & & & \\ \xi & \eta_{2} & \xi & 0 & \\ & & \ddots & & \\ & 0 & \xi & \eta_{2} & \xi \\ & & & \xi & \eta_{2} \end{pmatrix}\begin{pmatrix} u_{1,j-1}\\ u_{2,j-1}\\ \vdots\\ u_{m-2,j-1}\\ u_{m-1,j-1} \end{pmatrix}+\begin{pmatrix} -20\xi & \eta_{1} & & 0 & \\ \eta_{1} & -20\xi & \eta_{1} & & \\ & \ddots & \ddots & \ddots & \\ & & \eta_{1} & -20\xi & \eta_{1} \\ & 0 & & \eta_{1} & -20\xi \end{pmatrix}\begin{pmatrix} u_{1,j}\\ u_{2,j}\\ \vdots\\ u_{m-2,j}\\ u_{m-1,j} \end{pmatrix}+\begin{pmatrix} \eta_{2} & \xi & & & \\ \xi & \eta_{2} & \xi & 0 & \\ & & \ddots & & \\ & 0 & \xi & \eta_{2} & \xi \\ & & & \xi & \eta_{2} \end{pmatrix}\begin{pmatrix} u_{1,j+1}\\ u_{2,j+1}\\ \vdots\\ u_{m-2,j+1}\\ u_{m-1,j+1} \end{pmatrix}=\begin{pmatrix} -12\epsilon^{2}_{x}\epsilon^{2}_{y}f(x_{1},y_{j})-\xi (u_{0,j-1}+u_{0,j+1})-\eta_{1}u_{0,j}\\ -12\epsilon^{2}_{x}\epsilon^{2}_{y}f(x_{2},y_{j})\\ \vdots\\ -12\epsilon^{2}_{x}\epsilon^{2}_{y}f(x_{m-2},y_{j})\\ -12\epsilon^{2}_{x}\epsilon^{2}_{y}f(x_{m-1},y_{j})-\xi(u_{m,j-1}+u_{m,j+1})-\eta_{1}u_{m,j} \end{pmatrix}

其中,                          \xi=\beta+\gamma,\eta_{1}=10\beta-2\gamma,\eta_{2}=10\gamma-2\beta

记                                  \mathbf{u}_{j}=(u_{1,j},u_{2,j},\cdot\cdot\cdot,u_{m-1,j})^{T},0\leqslant j\leqslant n

则原数值格式可以简写为

 C(\mathbf{u}_{j-1}+\mathbf{u}_{j+1})+D\mathbf{u}_{j}=\mathbf{f}_{j},j=1,2,\cdot\cdot\cdot,n-1

其中,

C=\begin{pmatrix} \eta_{2} & \xi & & & \\ \xi & \eta_{2} & \xi & 0 & \\ & & \ddots & & \\ & 0 & \xi & \eta_{2} & \xi \\ & & & \xi & \eta_{2} \end{pmatrix}D=\begin{pmatrix} -20\xi & \eta_{1} & & 0 & \\ \eta_{1} & -20\xi & \eta_{1} & & \\ & \ddots & \ddots & \ddots & \\ & & \eta_{1} & -20\xi & \eta_{1} \\ & 0 & & \eta_{1} & -20\xi \end{pmatrix}

\mathbf{f}_{j}=\begin{pmatrix} -12\epsilon^{2}_{x}\epsilon^{2}_{y}f(x_{1},y_{j})-\xi (u_{0,j-1}+u_{0,j+1})-\eta_{1}u_{0,j}\\ -12\epsilon^{2}_{x}\epsilon^{2}_{y}f(x_{2},y_{j})\\ \vdots\\ -12\epsilon^{2}_{x}\epsilon^{2}_{y}f(x_{m-2},y_{j})\\ -12\epsilon^{2}_{x}\epsilon^{2}_{y}f(x_{m-1},y_{j})-\xi(u_{m,j-1}+u_{m,j+1})-\eta_{1}u_{m,j} \end{pmatrix}

        为解该方程组,将未知量\mathbf{u}_{j}按照下标拉长成一个列向量,并写出块矩阵形式,为

\begin{pmatrix} D & C & & & \\ C & D & C & & \\ & \ddots & \ddots & \ddots & \\ & & C & D & C\\ & & & C & D \end{pmatrix}\begin{pmatrix} \mathbf{u}_{1}\\ \mathbf{u}_{2}\\ \vdots\\ \mathbf{u}_{n-2}\\ \mathbf{u}_{n-1} \end{pmatrix}=\begin{pmatrix} \mathbf{f}_{1}-C\mathbf{u}_{0}\\ \mathbf{f}_{2}\\ \vdots\\ \mathbf{f}_{n-2}\\ \mathbf{f}_{n-1}-C\mathbf{u}_{n} \end{pmatrix}

上述方程组的高斯-赛德尔迭代公式为:

u^{(k+1)}_{i,j}=-\frac{1}{20\xi}[-12\epsilon^{2}_{x}\epsilon^{2}_{y}f(x_{i},y_{j})-\xi(u^{(k+1)}_{i-1,j-1}+u^{(k)}_{i-1,j+1}+u^{(k+1)}_{i+1,j-1}+u^{(k)}_{i+1,j+1})-\eta_{1}(u^{(k+1)}_{i-1,j}+u^{(k)}_{i+1,j})-\eta_{2}(u^{(k+1)}_{i,j-1}+u^{(k)}_{i,j+1})]\; (8)

其中,1\leqslant i \leqslant m-1, 1\leqslant j \leqslant n-1,加括号的上标k表示迭代次数。

三、算例实现

        采用九点紧差分格式求解椭圆型方程边值问题:

\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,\\ u(0,y)=sin(\pi y),u(2,y)=e^{2}sin(\pi y),0\leqslant y\leqslant 1,\\ u(x,0)=u(x,1)=0,0<x<2 \end{matrix}\right.

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

代码如下:(采用Gauss-Seidel迭代)


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,err,maxerr;double *x,*y,**u,**g,**temp,kexi,eta1,eta2;double leftboundary(double y);double rightboundary(double y);double bottomboundary(double x);double topboundary(double x);double f(double x, double y);double **Gij(double *x, double *y, int m, int n);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,m);dx=(xb-xa)/m;dy=(yb-ya)/n;beta=1.0/(dx*dx);gamma=1.0/(dy*dy);kexi=beta+gamma;eta1=10*beta-2*gamma;eta2=10*gamma-2*beta;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));temp=(double**)malloc(sizeof(double*)*(m+1));for(i=0;i<=m;i++){u[i]=(double*)malloc(sizeof(double)*(n+1));temp[i]=(double*)malloc(sizeof(double)*(n+1));}for(j=0;j<=n;j++){u[0][j]=leftboundary(y[j]);u[m][j]=rightboundary(y[j]);}for(i=1;i<m;i++){u[i][0]=bottomboundary(x[i]);u[i][n]=topboundary(x[i]);}for(i=1;i<m;i++){for(j=1;j<n;j++)u[i][j]=0.0;}g=Gij(x,y,m,n);for(i=0;i<=m;i++){for(j=0;j<=n;j++)temp[i][j]=u[i][j];}k=0;do{maxerr=0.0;for(i=1;i<m;i++){for(j=1;j<n;j++){temp[i][j]=(g[i][j]-kexi*(u[i-1][j-1]+temp[i-1][j+1]+u[i+1][j-1]+temp[i+1][j+1])-eta1*(u[i-1][j]+temp[i+1][j])-eta2*(u[i][j-1]+temp[i][j+1]))/(-20*kexi);err=temp[i][j]-u[i][j];if(err>maxerr)maxerr=err;u[i][j]=temp[i][j];}}k=k+1;}while(maxerr>0.5*1e-10);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(temp[i]);}free(u);free(temp);free(x);free(y);return 0;
}double leftboundary(double y)
{return sin(pi*y);
}
double rightboundary(double y)
{return exp(1.0)*exp(1.0)*sin(pi*y);
}
double bottomboundary(double x)
{return 0.0;
}
double topboundary(double x)
{return 0.0;
}
double exact(double x, double y)
{return exp(x)*sin(pi*y);
}
double f(double x, double y)
{return (pi*pi-1)*exp(x)*sin(pi*y);
}
double **Gij(double *x, double *y, int m, int n)
{int i,j;double temp1,temp2,temp3,**ans;ans=(double**)malloc(sizeof(double*)*(m+1));for(i=0;i<=m;i++)ans[i]=(double*)malloc(sizeof(double)*(n+1));for(i=1;i<m;i++){for(j=1;j<n;j++){temp1=f(x[i-1],y[j-1])+10*f(x[i],y[j-1])+f(x[i+1],y[j-1]);temp2=f(x[i-1],y[j])+10*f(x[i],y[j])+f(x[i+1],y[j]);temp3=f(x[i-1],y[j+1])+10*f(x[i],y[j+1])+f(x[i+1],y[j+1]);ans[i][j]=-(temp1+temp3+10*temp2)/12.0;}}return ans;
}

 相同的问题,与五点菱形差分法计算结果进行了对比,结果如下:

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

+++++++++++++九点紧差分格式++++++++++++++
m=32, n=16.
k=747.
(0.50,0.25), y=1.165829, err=6.7140e-06.
(1.00,0.25), y=1.922127, err=1.1487e-05.
(1.50,0.25), y=3.169047, err=1.4319e-05.
(0.50,0.50), y=1.648731, err=9.4951e-06.
(1.00,0.50), y=2.718298, err=1.6245e-05.
(1.50,0.50), y=4.481709, err=2.0250e-05.

+++++++++++++五点菱形差分格式++++++++++++
m=32,n=16.
k=887
(0.50,0.25), y=1.169343, err=3.5207e-03.
(1.00,0.25), y=1.928138, err=6.0228e-03.
(1.50,0.25), y=3.176531, err=7.4986e-03.
(0.50,0.50), y=1.653700, err=4.9790e-03.
(1.00,0.50), y=2.726799, err=8.5175e-03.
(1.50,0.50), y=4.492294, err=1.0605e-02.

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

+++++++++++++九点紧差分格式++++++++++++++
m=64, n=32.
k=2790.
(0.50,0.25), y=1.165822, err=4.1552e-07.
(1.00,0.25), y=1.922116, err=7.1227e-07.
(1.50,0.25), y=3.169034, err=8.9066e-07.
(0.50,0.50), y=1.648722, err=5.8773e-07.
(1.00,0.50), y=2.718283, err=1.0074e-06.
(1.50,0.50), y=4.481690, err=1.2597e-06.

 

+++++++++++++五点菱形差分格式++++++++++++
m=64,n=32.
k=3315
(0.50,0.25), y=1.166702, err=8.7958e-04.
(1.00,0.25), y=1.923620, err=1.5048e-03.
(1.50,0.25), y=3.170908, err=1.8751e-03.
(0.50,0.50), y=1.649965, err=1.2439e-03.
(1.00,0.50), y=2.720410, err=2.1281e-03.
(1.50,0.50), y=4.484341, err=2.6518e-03.

四、结论

        从计算结果可知,对于九点紧差分格式来说,当步长减半时,误差减小为约1/16,可见其数值格式是四阶收敛的。同时,与五点菱形差分格式相比,九点紧差分法迭代步数更少,同时计算精度更高!

这篇关于偏微分方程算法之九点紧差分法的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



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

相关文章

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

本文以商业化应用推荐为例,告诉我们不懂推荐算法的产品,也能从产品侧出发, 设计出一款不错的推荐系统。 相信很多新手产品,看到算法二字,多是懵圈的。 什么排序算法、最短路径等都是相对传统的算法(注:传统是指科班出身的产品都会接触过)。但对于推荐算法,多数产品对着网上搜到的资源,都会无从下手。特别当某些推荐算法 和 “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%免费

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份,对于每一份内部的数进行翻转(逆序),每次操作完后输出操作后新序列的逆序对数。 图一:  划分子问题。 图二: 分而治之,=>  合并 。 图三: 回溯:

最大公因数:欧几里得算法

简述         求两个数字 m和n 的最大公因数,假设r是m%n的余数,只要n不等于0,就一直执行 m=n,n=r 举例 以18和12为例 m n r18 % 12 = 612 % 6 = 06 0所以最大公因数为:6 代码实现 #include<iostream>using namespace std;/