本文主要是介绍偏微分方程算法之二阶双曲型方程交替方向隐格式,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!
目录
一、研究目标
二、理论推导
2.1 显差分格式
2.2 交替方向隐格式
三、算例实现
四、结论
一、研究目标
本节的研究对象为如下的双曲型偏微分方程初边值问题:
为了保证连续性,公式(1)中相关函数满足:
二、理论推导
2.1 显差分格式
第一步:网格剖分。在三维长方体空间(二维平面加一维时间轴)进行网格剖分,将区域[0,a]等分m份,区域[0,b]等分n份,区域[0,T]等分l份,即
得到网格节点坐标为。我们利用数值方法获得精确解在网格节点处的近似值,即数值解。
第二步:弱化原方程,仅在离散点处成立。即
第三步:差商代替微商处理偏导数。方程组二阶偏导数使用中心差商来近似,即
对于初始条件中的一阶偏导数,也采用二阶中心差商,即
将差商公式带入公式(2),用数值解代替精确解并忽略高阶项,可得对应的数值格式为
为消去虚拟边界项,让公式(3)中第1式对k=0成立,可得
再结合公式(3)中第3式中越界的,消去虚拟边界项,可得
记,整理上式可得
由公式(4)知,公式(3)即为三层显格式:
该格式的局部截断误差为,考虑到其稳定性限制,将其改为隐格式。
2.2 交替方向隐格式
由于
由上面两式可将原来在点处的连续方程
改写为
再用中心差商处理二阶偏导数,将数值解代替精确解并忽略高阶项,可得数值隐格式为:
公式(5)在求解上存在困难,故参照二维抛物型偏微分方程初边值问题的交替方向隐格式设计思想,先将公式(5)写成算子的形式,可改写为:
也就是
即
对上式左端添加一个辅助项,即
对公式(6)进行分解,可以写为
通过引入中间变量,得
上式是一个x方向的三对角线性方程组,计算第1式中的时需要用到和,这可以直接从第2式中解得,即
我们可以直接将公式(7)、(8)写成差分形式,即
其中,
假设第k和k-1个时间层的信息已知,这样在公式(9)的第1式中,先固定某个j,这时公式(9)的第1式就是一个m-1阶线性方程组,其系数矩阵是三对角矩阵,可以用追赶法求解。每解一个方程组就可以得到一个在[0,a]x[0,b]矩形区域内行网格点(即行节点)上的中间量V的信息,所以要得到矩形区域内所有的网格点信息,需要求解n-1个这样的线性系统。同样,公式(9)的第2式是一个y方向的三对角线性方程组时,求解时要先固定某个i,这时公式(9)的第2式是一个n-1阶线性方程组,系数矩阵是三对角矩阵,可以用追赶法求解,每解一个方程组得到第k+1个时间层上[0,a]x[0,b]矩形区域内列网格点(即)上的数值解的信息,所以要得到第k+1时间层上所有网格点的信息,需要求解m-1个这样的线性系统。这样就实现了x方向和y方向交替求解的隐格式。
三、算例实现
用交替方向隐格式求解下列二维双曲型方程初边值问题:
已知精确解为。分别取步长为和,给出在节点处的数值解和误差。
代码如下:
#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);
}
当时,计算结果如下:
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
当时,计算结果如下:
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,可见该数值格式是二阶的。
这篇关于偏微分方程算法之二阶双曲型方程交替方向隐格式的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!