Пример. Уравнение теплопроводности

Или как работать с изменением по времени?

\[\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')