考虑如下两点边值问题
\[\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()