2024 XTU PDE 12题 两点边值问题有限体积格式

考虑如下两点边值问题
\[\begin{cases}-u^{''}+10(2x-1)u'+20u=0,\quad0<x<1,\\
u(0)=u(1)=1.\end{cases}\]
真解取为
\[u=e^{-10x(1-x)}.\]
写出求解上述问题的有限体积格式。
编程实现上述有限体积法，并计算节点处的最大误差，考查其收敛阶。

解：(a) 有限体积格式：
\[\left\{
\begin{aligned}
&-\left[\frac{u_{i+1}-u_i}{h_{i+1}}-
\frac{u_i-u_{i-1}}{h_i}\right]+\frac12c_i\left(u_{i+1}-u_{i-1}\right)+\frac{h_i+h_{i+1}}2d_iu_i=0,\quad i=1,2,\cdots,N-1.\\
&u_0=u_N=1.
\end{aligned}\right.
\]
其中
\[ c_i=\frac2{h_i+h_{i+1}}\int_{x_{i-1/2}}^{x_{i+1/2}}10(2x-1)dx,\quad d_i=\frac2{h_i+h_{i+1}}\int_{x_{i-1/2}}^{x_{i+1/2}}20dx.\]

import numpy as np
from scipy.sparse import coo_matrix, csc_matrix, \
    csr_matrix, spdiags
from scipy.sparse.linalg import spsolve
import matplotlib.pyplot as plt


class PoissonData1d:
    def __init__(self, L=0, R=1):
        self.L = L
        self.R = R

    def domain(self):
        return (self.L, self.R)

    def init_mesh(self, NS):
        x = np.linspace(self.L, self.R, num=NS + 1)
        return x

    def solution(self, x):
        """
        :param x:
        :return: 真解
        """
        return np.exp(-10 * x * (1 - x))

    def source(self, x):
        """
        :param x:
        :return: 右端向量
        """
        return np.zeros(len(x))


def fvm(data, NS):
    x = data.init_mesh(NS=NS)
    NN = len(x)
    h = x[1] - x[0]
    c1 = np.zeros(NN - 2)
    for i in range(NN - 2):
        c1[i] = -1 / h + 10 * x[i+1] - 5
    val = np.zeros((3, NN - 2))
    c2 = 2 / h + 20 * h
    c3 = np.zeros(NN - 2)
    for i in range(NN - 2):
        c3[i] = -1 / h - 10 * x[i+1] + 5
    val[1, :] = c2
    val[0, :] = c1

    val[2, :] = c3

    diags = np.array([-1, 0, 1])
    A = spdiags(val, diags, NN - 2, NN - 2).tocsr()
    A = A.T
    # print(A)
    F = np.zeros(len(x[1: NN - 1]))
    F[0] = -c3[0]
    F[-1] = -c1[-1]
    print(F)
    u = np.zeros((NN,), dtype=np.float_)
    u[[0, -1]] = data.solution(x[[0, -1]])
    # print(A)
    u[1: -1] = spsolve(A, F)
    return x, u


def error(x, u, uh):
    e = u - uh
    emax = np.max(np.abs(e))
    return emax


L = 0
R = 1
NS = 4
maxit = 5
data = PoissonData1d(L=L, R=R)
e = np.zeros((maxit, 1), dtype=np.float_)
fig = plt.figure()
axes = fig.gca()
for i in range(maxit):
    x, uh = fvm(data, NS)
    u = data.solution(x)
    e[i, 0] = error(x, u, uh)
    axes.plot(x, uh, linestyle='-', marker='o', markersize=4, label='N=%d' % NS)
    axes.set_title(" The solution plot ")
    axes.set_xlabel("x")
    axes.set_ylabel("u")
    if i < maxit - 1:
        NS *= 2

axes.plot(x, u, label='exact')
plt.legend(loc='upper right')
print(" error :\n", e)
plt.show()