Пример. Уравнение теплопроводности#
Или как работать с изменением по времени?
\[\begin{split}
\left\{\begin{array}{lll}
\frac{\partial u}{\partial t}=\alpha \frac{\partial^2 u}{\partial x^2} & \text { in } \Omega=[0,1] & t>0 \\
u=u_0 & \text { on } \partial \Omega & t>0 \\
u=I & \text { in } \Omega & t=0
\end{array}\right.
\end{split}\]
Дискретизуем уравнение по времени, а потом решаем с помощью конечных элементов.
\[
\frac{u^{[t+1]}-u^{[t]}}{\Delta t}=\alpha \frac{\partial^2}{\partial x^2} u^{[t+1]}
\]
Приведем его к слабой форме.
\[\begin{split}
\begin{aligned}
& \int_{\Omega}\left(u^{[t+1]}-u^{[t]}\right) v=\alpha \Delta t \int_{\Omega} \frac{\partial^2}{\partial x^2} u^{[t+1]} v \\
& \int_{\Omega}\left(u^{[t+1]}-u^{[t]}\right) v=-\alpha \Delta t \int_{\Omega} \frac{\partial}{\partial x} u^{[t+1]} \frac{\partial}{\partial x} v . \\
& \int_{\Omega} u^{[t+1]} v-\alpha \Delta t \nabla u^{[t+1]} \nabla v d x=\int_{\Omega} u^{[t]} v d x
\end{aligned}
\end{split}\]
Уравнение имеет вид:
Билинейная форма \(\left(u^{[t+1]}, v\right)=\) Линейная форма \((v)\) (Мы знаем \(u^{[t]}\) с предыдущего шаг).
N_elements = 20
mesh = fe.UnitIntervalMesh(N_elements)
lagr_pol_ford = fe.FunctionSpace(
mesh,
'Lagrange',
1
)
u_on_b_l = fe.Constant(0.0)
u_on_b_r = fe.Constant(-2.0)
def near(x, const, tol):
return abs(x - const) < tol
def boundary_L(x, on_boundary):
tol = 1e-14
return on_boundary and near(x[0], 0, tol)
def boundary_R(x, on_boundary):
tol = 1e-14
return on_boundary and near(x[0], 1, tol)
bs_l = fe.DirichletBC(
lagr_pol_ford,
u_on_b_l,
boundary_L,
)
bs_r = fe.DirichletBC(
lagr_pol_ford,
u_on_b_r,
boundary_R,
)
bs = [bs_l, bs_r]
initial_condition = fe.Expression(
'sin(3.1415 * x[0]) * x[0] * x[0] * exp(2*x[0]) - 2*x[0]',
degree = 1
)
u_old = fe.interpolate(
initial_condition,
lagr_pol_ford
)
dt = 0.1
heat_source = fe.Constant(0.0)
u_trial = fe.TrialFunction(lagr_pol_ford)
v_test = fe.TestFunction(lagr_pol_ford)
weak_form_residuum = (
u_trial * v_test * fe.dx
+
dt * fe.dot(
fe.grad(u_trial),
fe.grad(v_test),
) * fe.dx
-
(
u_old * v_test * fe.dx
+
dt * heat_source * v_test * fe.dx
)
)
weak_form_lhs = fe.lhs(weak_form_residuum)
weak_form_rhs = fe.rhs(weak_form_residuum)
plt.figure()
fe.plot(u_old, label='time=0.0')
n_time_steps = 5
time_current = 0
u_solution = fe.Function(lagr_pol_ford)
for i in range(n_time_steps):
time_current += dt
fe.solve(
weak_form_lhs == weak_form_rhs,
u_solution,
bs,
)
u_old.assign(u_solution)
fe.plot(u_solution, label = 'time = {:.2f}'.format(time_current))
plt.legend()
plt.xlabel('x')
plt.ylabel('u')
Text(0, 0.5, 'u')