Phương trình nhiệt 3D Python

Là phần mở rộng đầu tiên của bài toán Poisson từ chương trước, chúng ta xét phương trình nhiệt phụ thuộc thời gian, hay phương trình khuếch tán phụ thuộc thời gian. Đây là phần mở rộng tự nhiên của phương trình Poisson mô tả sự phân bố nhiệt cố định trong cơ thể đối với một vấn đề phụ thuộc vào thời gian

Chúng ta sẽ thấy rằng bằng cách rời rạc hóa thời gian thành các khoảng thời gian nhỏ và áp dụng các phương pháp bước thời gian tiêu chuẩn, chúng ta có thể giải phương trình nhiệt bằng cách giải một chuỗi các bài toán biến phân, giống như bài toán chúng ta đã gặp đối với phương trình Poisson

vấn đề PDE

Bài toán mô hình của chúng tôi đối với các PDE phụ thuộc vào thời gian có dạng $$ \begin{alignat}{2} {\partial u\over\partial t} &= \nabla^2 u + f \quad &&\hbox{in }\Omega\times . 1}\\ u &= \ub &&\hbox{on } \partial \Omega\times[0, T], \tag{3. 2}\\ u &= \I &\m hộp{at } t=0\to \tag{3. 3} \end{alignat} $$ Ở đây, \[ u \] thay đổi theo không gian và thời gian, e. g. , \[ u=u[x,y,t] \] nếu miền không gian \[ \Omega \] là hai chiều. Hàm nguồn \[ f \] và các giá trị biên \[ \ub \] cũng có thể thay đổi theo không gian và thời gian. Điều kiện ban đầu \[ \uI \] chỉ là một hàm của không gian

Một cách tiếp cận đơn giản để giải các PDE phụ thuộc vào thời gian bằng phương pháp phần tử hữu hạn là trước tiên rời rạc hóa đạo hàm thời gian bằng một xấp xỉ sai phân hữu hạn, tạo ra một chuỗi các bài toán dừng, sau đó biến từng bài toán dừng thành một công thức biến phân

Đặt chỉ số trên \[ n \] biểu thị một đại lượng tại thời điểm \[ t_n \], trong đó \[ n \] là một số nguyên đếm các mức thời gian. Ví dụ: \[ u^n \] có nghĩa là \[ u \] ở mức thời gian \[ n \]. Trước tiên, một sự khác biệt hữu hạn rời rạc theo thời gian bao gồm lấy mẫu PDE ở một mức thời gian nào đó, giả sử \[ t_{n+1} \]. $$ \begin{equation} \left[{\partial u \over\partial t}\right]^{n+1} = \nabla^2 u^{n+1} + f^{n+1}\ . 4} \end{equation} $$ Đạo hàm thời gian có thể xấp xỉ bằng thương số. Vì lý do đơn giản và ổn định, ta chọn sai phân lùi đơn giản. $$ \begin{equation} \left[{\partial u\over\partial t}\right]^{n+1}\approx {{u^{n+1} - u^n}\over{\dt . 5} \end{equation} $$ trong đó \[ \dt \] là tham số rời rạc về thời gian. Chèn vào kết quả là $$ \begin{equation} {{u^{n+1} - u^n}\over{\dt}} = \nabla^2 u^{n+1} + f^{n+1 . 6} \end{equation} $$ Đây là phiên bản rời rạc theo thời gian của phương trình nhiệt , cái gọi là Euler ngược hoặc rời rạc Euler ẩn

Chúng tôi có thể sắp xếp lại để phía bên tay trái chứa các thuật ngữ chưa biết \[ u^{n+1} \] và phía bên tay phải chỉ chứa các thuật ngữ được tính toán. Kết quả là một chuỗi các bài toán không gian [cố định] cho \[ u^{n+1} \], giả sử \[ u^n \] đã biết từ bước thời gian trước đó. $$ \begin{align} u^0 &= \uI, \tag{3. 7}\\ u^{n+1} - {\dt}\nabla^2 u^{n+1} &= u^n + {\dt} f^{n+1},\quad n=0 . 8} \end{align} $$ Với \[ \uI \], chúng ta có thể giải ra \[ u^0 \], \[ u^1 \], \[ u^2 \], v.v.

Một thay thế cho , có thể thuận tiện trong việc triển khai, là thu thập tất cả các số hạng ở một bên của dấu đẳng thức. $$ \begin{equation} u^{n+1} - {\dt}\nabla^2 u^{n+1} - u^{n} - {\dt} f^{n+1} = 0 . 9} \end{phương trình} $$

Chúng tôi sử dụng một phương pháp phần tử hữu hạn để giải quyết và một trong hai phương trình hoặc. Điều này đòi hỏi phải chuyển các phương trình về dạng yếu. Như thường lệ, chúng tôi nhân với hàm kiểm tra \[ v\in \hat V \] và tích phân đạo hàm cấp hai theo từng phần. Giới thiệu ký hiệu \[ u \] cho \[ u^{n+1} \] [điều này là tự nhiên trong chương trình], dạng yếu phát sinh từ công thức có thể được viết một cách thuận tiện trong ký hiệu tiêu chuẩn. $$ a[u,v]=L_{n+1}[v],$$ trong đó $$ \begin{align} a[u,v] &= \int_\Omega\left[uv + {\dt} . 10}\\ L_{n+1}[v] &= \int_\Omega \left[u^n + {\dt} f^{n+1}\right]v \dx\tp \tag{3. 11} \end{align} $$ Dạng thay thế có công thức trừu tượng $$ F_{n+1}[u;v] = 0,$$ trong đó $$ \begin{equation} F_{n+1}[u . 12} \end{phương trình} $$

Ngoài bài toán biến phân cần giải theo từng bước thời gian, ta còn cần xấp xỉ điều kiện ban đầu. Phương trình này cũng có thể được chuyển thành một bài toán biến phân. $$ a_0[u,v]=L_0[v],$$ với $$ \begin{align} a_0[u,v] &= \int_\Omega uv \dx, \tag{3. 13}\\ L_0[v] &= \int_\Omega \uI v \dx\tp \tag{3. 14} \end{align} $$ Khi giải bài toán biến phân này, \[ u^0 \] trở thành phép chiếu \[ L^2 \] của giá trị ban đầu \[ \uI \] đã cho vào không gian phần tử hữu hạn. Cách khác là xây dựng \[ u^0 \] bằng cách chỉ nội suy giá trị ban đầu \[ \uI \]; . Chúng tôi gọi hai chiến lược này là tính toán điều kiện ban đầu bằng phép chiếu hoặc phép nội suy. Cả hai hoạt động đều dễ dàng tính toán trong FEniCS thông qua một câu lệnh duy nhất, sử dụng hàm

def boundary[x, on_boundary]:
    return on_boundary

bc = DirichletBC[V, u_D, boundary]

0 hoặc

def boundary[x, on_boundary]:
    return on_boundary

bc = DirichletBC[V, u_D, boundary]

1. Lựa chọn phổ biến nhất là

def boundary[x, on_boundary]:
    return on_boundary

bc = DirichletBC[V, u_D, boundary]

0, tính toán gần đúng với \[ \uI \], nhưng trong một số ứng dụng mà chúng tôi muốn xác minh mã bằng cách sao chép các giải pháp chính xác, người ta phải sử dụng

def boundary[x, on_boundary]:
    return on_boundary

bc = DirichletBC[V, u_D, boundary]

1 [và chúng tôi sử dụng bài toán kiểm tra như vậy tại đây. ]

Tóm lại, do đó chúng ta cần giải chuỗi các bài toán biến phân sau đây để tính nghiệm phần tử hữu hạn cho phương trình nhiệt. tìm \[ u^0\in V \] sao cho \[ a_0[u^0,v]=L_0[v] \] đúng với mọi \[ v\in\hat V \], rồi tìm \[ u

Chương trình của chúng tôi cần triển khai bước thời gian theo cách thủ công nhưng có thể dựa vào FEniCS để dễ dàng tính toán \[ a_0 \], \[ L_0 \], \[ a \] và \[ L \] [hoặc \[ F_{n

Đề kiểm tra 1. Một giải pháp phân tích đã biết

Cũng giống như bài toán Poisson ở chương trước, chúng ta xây dựng một bài toán kiểm tra giúp dễ dàng xác định xem các phép tính có đúng không. Vì chúng tôi biết rằng sơ đồ bước thời gian bậc nhất của chúng tôi là chính xác cho các hàm tuyến tính, chúng tôi tạo một bài toán kiểm tra có biến thiên tuyến tính theo thời gian. Chúng tôi kết hợp điều này với một biến thể bậc hai trong không gian. Do đó, chúng ta lấy $$ \begin{equation} u = 1 + x^2 + \alpha y^2 + \beta t, \tag{3. 15} \end{equation} $$ tạo ra một hàm có giá trị được tính toán tại các nút sẽ chính xác, bất kể kích thước của các phần tử và \[ \dt \], miễn là lưới được phân chia đồng đều. Bằng cách chèn vào phương trình nhiệt , chúng ta thấy rằng vế phải \[ f \] phải được cho bởi \[ f[x,y,t]=\beta - 2 - 2\alpha \]. Giá trị biên là \[ \ub[x, y, t] = 1 + x^2 + \alpha y^2 + \beta t \] và giá trị ban đầu là \[ \uI[x, y] = 1 +

triển khai FEniCS

Một vấn đề lập trình mới là cách xử lý các hàm thay đổi theo không gian và thời gian, chẳng hạn như điều kiện biên \[ \ub[x, y, t] = 1 + x^2 + \alpha y^2 + \beta t \ . Một giải pháp tự nhiên là sử dụng FEniCS

def boundary[x, on_boundary]:
    return on_boundary

bc = DirichletBC[V, u_D, boundary]

4 với thời gian \[ t \] làm tham số, ngoài các tham số \[ \alpha \] và \[ \beta \]

alpha = 3; beta = 1.2
u_D = Expression['1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t',
                 degree=2, alpha=alpha, beta=beta, t=0]

def boundary[x, on_boundary]:
    return on_boundary

bc = DirichletBC[V, u_D, boundary]

4 này sử dụng các thành phần của

def boundary[x, on_boundary]:
    return on_boundary

bc = DirichletBC[V, u_D, boundary]

6 làm biến độc lập, trong khi

def boundary[x, on_boundary]:
    return on_boundary

bc = DirichletBC[V, u_D, boundary]

7,

def boundary[x, on_boundary]:
    return on_boundary

bc = DirichletBC[V, u_D, boundary]

8 và

def boundary[x, on_boundary]:
    return on_boundary

bc = DirichletBC[V, u_D, boundary]

9 là các tham số. Thời gian

def boundary[x, on_boundary]:
    return on_boundary

bc = DirichletBC[V, u_D, boundary]

9 sau này có thể được cập nhật bởi

Các điều kiện biên thiết yếu, dọc theo toàn bộ đường biên trong trường hợp này, được triển khai giống như cách chúng ta đã triển khai các điều kiện biên cho bài toán Poisson trước đây

def boundary[x, on_boundary]:
    return on_boundary

bc = DirichletBC[V, u_D, boundary]

Chúng ta sẽ sử dụng biến

u_n = project[u_D, V]
# or
u_n = interpolate[u_D, V]

1 cho \[ u^{n+1} \] chưa biết ở bước thời gian mới và biến

u_n = project[u_D, V]
# or
u_n = interpolate[u_D, V]

2 cho \[ u^n \] ở bước thời gian trước đó. Giá trị ban đầu của

u_n = project[u_D, V]
# or
u_n = interpolate[u_D, V]

2 có thể được tính bằng phép chiếu hoặc phép nội suy của \[ \uI \]. Vì chúng tôi đặt

u_n = project[u_D, V]
# or
u_n = interpolate[u_D, V]

4 cho giá trị biên

u_n = project[u_D, V]
# or
u_n = interpolate[u_D, V]

5, nên chúng tôi có thể sử dụng

u_n = project[u_D, V]
# or
u_n = interpolate[u_D, V]

5 để chỉ định điều kiện ban đầu

u_n = project[u_D, V]
# or
u_n = interpolate[u_D, V]

Chiếu so với nội suy điều kiện ban đầu. Để thực sự khôi phục giải pháp chính xác cho độ chính xác của máy, điều quan trọng là phải tính toán điều kiện ban đầu rời rạc bằng cách nội suy \[ \uI \]. Điều này đảm bảo rằng bậc tự do là chính xác [với độ chính xác của máy] tại \[ t=0 \]. Phép chiếu dẫn đến các giá trị gần đúng tại các nút

Chúng tôi có thể xác định \[ a \] hoặc \[ L \] theo các công thức trên hoặc chúng tôi có thể chỉ xác định \[ F \] và yêu cầu FEniCS tìm ra các thuật ngữ nào sẽ chuyển sang dạng song tuyến tính \[ a \] . Cái sau thuận tiện, đặc biệt là trong các vấn đề phức tạp hơn, vì vậy chúng tôi minh họa cách xây dựng \[ a \] và \[ L \]

u = TrialFunction[V]
v = TestFunction[V]
f = Constant[beta - 2 - 2*alpha]

F = u*v*dx + dt*dot[grad[u], grad[v]]*dx - [u_n + dt*f]*v*dx
a, L = lhs[F], rhs[F]

Cuối cùng, chúng tôi thực hiện bước thời gian trong một vòng lặp

u = Function[V]
t = 0
for n in range[num_steps]:

    # Update current time
    t += dt
    u_D.t = t

    # Solve variational problem
    solve[a == L, u, bc]

    # Update previous solution
    u_n.assign[u]

Ở bước cuối cùng của vòng lặp bước thời gian, chúng ta gán các giá trị của biến

u_n = project[u_D, V]
# or
u_n = interpolate[u_D, V]

1 [nghiệm pháp được tính toán mới] cho biến

u_n = project[u_D, V]
# or
u_n = interpolate[u_D, V]

2 chứa các giá trị ở bước thời gian trước đó. Điều này phải được thực hiện bằng cách sử dụng hàm thành viên

u_n = project[u_D, V]
# or
u_n = interpolate[u_D, V]

9. Thay vào đó, nếu chúng tôi cố gắng thực hiện

u = TrialFunction[V]
v = TestFunction[V]
f = Constant[beta - 2 - 2*alpha]

F = u*v*dx + dt*dot[grad[u], grad[v]]*dx - [u_n + dt*f]*v*dx
a, L = lhs[F], rhs[F]

0, chúng tôi sẽ đặt biến

u_n = project[u_D, V]
# or
u_n = interpolate[u_D, V]

2 thành biến giống như biến

u_n = project[u_D, V]
# or
u_n = interpolate[u_D, V]

1, đây không phải là điều chúng tôi muốn. [Chúng tôi cần hai biến, một cho các giá trị ở bước thời gian trước và một cho các giá trị ở bước thời gian hiện tại. ]

Hãy nhớ cập nhật các đối tượng biểu thức với thời gian hiện tại. Bên trong vòng lặp thời gian, hãy quan sát rằng _____20_______3 phải được cập nhật trước câu lệnh ____20_______4 để thực thi tính toán các điều kiện Dirichlet tại bước thời gian hiện tại. Một điều kiện Dirichlet được xác định dưới dạng một

def boundary[x, on_boundary]:
    return on_boundary

bc = DirichletBC[V, u_D, boundary]

4 tra cứu và áp dụng giá trị của một tham số chẳng hạn như

def boundary[x, on_boundary]:
    return on_boundary

bc = DirichletBC[V, u_D, boundary]

9 khi nó được đánh giá và áp dụng cho hệ thống tuyến tính

Vòng lặp bước thời gian ở trên không chứa bất kỳ so sánh nào về các giải pháp số và giải pháp chính xác mà chúng tôi phải đưa vào để xác minh việc triển khai. Đối với phương trình Poisson trong phần , chúng tôi tính toán sự khác biệt giữa mảng các giá trị nốt cho

u_n = project[u_D, V]
# or
u_n = interpolate[u_D, V]

1 và mảng các giá trị nốt cho nghiệm chính xác được nội suy. Điều này có thể được thực hiện như sau

u_e = interpolate[u_D, V]
error = np.abs[u_e.vector[].array[] - u.vector[].array[]].max[]
print['t = %.2f: error = %.3g' % [t, error]]

Đối với ví dụ Poisson, chúng tôi đã sử dụng hàm

u = TrialFunction[V]
v = TestFunction[V]
f = Constant[beta - 2 - 2*alpha]

F = u*v*dx + dt*dot[grad[u], grad[v]]*dx - [u_n + dt*f]*v*dx
a, L = lhs[F], rhs[F]

8 để trích xuất các giá trị hàm tại các đỉnh. Ở đây chúng tôi minh họa một phương pháp khác để trích xuất các giá trị đỉnh, bằng cách gọi hàm

u = TrialFunction[V]
v = TestFunction[V]
f = Constant[beta - 2 - 2*alpha]

F = u*v*dx + dt*dot[grad[u], grad[v]]*dx - [u_n + dt*f]*v*dx
a, L = lhs[F], rhs[F]

9, trả về vectơ bậc tự do. Đối với không gian hàm \[ \mathsf{P}_1 \], vectơ bậc tự do này sẽ bằng với mảng các giá trị đỉnh thu được bằng cách gọi

u = TrialFunction[V]
v = TestFunction[V]
f = Constant[beta - 2 - 2*alpha]

F = u*v*dx + dt*dot[grad[u], grad[v]]*dx - [u_n + dt*f]*v*dx
a, L = lhs[F], rhs[F]

8, mặc dù có thể theo thứ tự khác

Chương trình hoàn chỉnh để giải phương trình nhiệt diễn ra như sau

from fenics import *
import numpy as np

T = 2.0            # final time
num_steps = 10     # number of time steps
dt = T / num_steps # time step size
alpha = 3          # parameter alpha
beta = 1.2         # parameter beta

# Create mesh and define function space
nx = ny = 8
mesh = UnitSquareMesh[nx, ny]
V = FunctionSpace[mesh, 'P', 1]

# Define boundary condition
u_D = Expression['1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t',
                 degree=2, alpha=alpha, beta=beta, t=0]

def boundary[x, on_boundary]:
    return on_boundary

bc = DirichletBC[V, u_D, boundary]

# Define initial value
u_n = interpolate[u_D, V]
#u_n = project[u_D, V]

# Define variational problem
u = TrialFunction[V]
v = TestFunction[V]
f = Constant[beta - 2 - 2*alpha]

F = u*v*dx + dt*dot[grad[u], grad[v]]*dx - [u_n + dt*f]*v*dx
a, L = lhs[F], rhs[F]

# Time-stepping
u = Function[V]
t = 0
for n in range[num_steps]:

    # Update current time
    t += dt
    u_D.t = t

    # Compute solution
    solve[a == L, u, bc]

    # Plot solution
    plot[u]

    # Compute error at vertices
    u_e = interpolate[u_D, V]
    error = np.abs[u_e.vector[].array[] - u.vector[].array[]].max[]
    print['t = %.2f: error = %.3g' % [t, error]]

    # Update previous solution
    u_n.assign[u]

# Hold plot
interactive[]

Chương trình ví dụ này có thể được tìm thấy trong tệp ft03_heat. py

Đề kiểm tra 2. Sự khuếch tán của hàm Gaussian

Bây giờ chúng ta hãy giải một bài toán kiểm tra thú vị hơn, đó là sự khuếch tán của một ngọn đồi Gaussian. Chúng tôi lấy giá trị ban đầu là $$ \uI[x,y]= e^{-ax^2 - ay^2}$$ cho \[ a = 5 \] trên miền \[ [-2,2] . Đối với vấn đề này, chúng tôi sẽ sử dụng các điều kiện biên Dirichlet thuần nhất [\[ \ub = 0 \]]

triển khai FEniCS

Đó là những thay đổi bắt buộc đối với chương trình trước đây của chúng tôi? . Miền mới có thể được tạo dễ dàng trong FEniCS bằng cách sử dụng

u = Function[V]
t = 0
for n in range[num_steps]:

    # Update current time
    t += dt
    u_D.t = t

    # Solve variational problem
    solve[a == L, u, bc]

    # Update previous solution
    u_n.assign[u]

1

nx = ny = 30
mesh = RectangleMesh[Point[-2, -2], Point[2, 2], nx, ny]

Lưu ý rằng chúng tôi đã sử dụng độ phân giải cao hơn nhiều so với trước đây để giải quyết tốt hơn các tính năng của giải pháp. Ta cũng cần xác định lại điều kiện ban đầu và điều kiện biên. Cả hai đều dễ dàng thay đổi bằng cách xác định một

def boundary[x, on_boundary]:
    return on_boundary

bc = DirichletBC[V, u_D, boundary]

4 mới và bằng cách đặt \[ u = 0 \] trên ranh giới

Để có thể trực quan hóa giải pháp trong một chương trình bên ngoài như ParaView, chúng tôi sẽ lưu giải pháp vào một tệp ở định dạng VTK trong mỗi bước thời gian. Chúng tôi thực hiện việc này bằng cách tạo một

u = Function[V]
t = 0
for n in range[num_steps]:

    # Update current time
    t += dt
    u_D.t = t

    # Solve variational problem
    solve[a == L, u, bc]

    # Update previous solution
    u_n.assign[u]

3 đầu tiên với hậu tố

u = Function[V]
t = 0
for n in range[num_steps]:

    # Update current time
    t += dt
    u_D.t = t

    # Solve variational problem
    solve[a == L, u, bc]

    # Update previous solution
    u_n.assign[u]

4

vtkfile = File['heat_gaussian/solution.pvd']

Trong vòng lặp thời gian, sau đó chúng tôi có thể nối các giá trị giải pháp vào tệp này

Dòng này được gọi trong mỗi bước thời gian, dẫn đến việc tạo một tệp mới có hậu tố

u = Function[V]
t = 0
for n in range[num_steps]:

    # Update current time
    t += dt
    u_D.t = t

    # Solve variational problem
    solve[a == L, u, bc]

    # Update previous solution
    u_n.assign[u]

5 chứa tất cả dữ liệu cho bước thời gian [giá trị lưới và đỉnh]. Tệp

u = Function[V]
t = 0
for n in range[num_steps]:

    # Update current time
    t += dt
    u_D.t = t

    # Solve variational problem
    solve[a == L, u, bc]

    # Update previous solution
    u_n.assign[u]

6 sẽ chứa các giá trị thời gian và tham chiếu đến tệp

u = Function[V]
t = 0
for n in range[num_steps]:

    # Update current time
    t += dt
    u_D.t = t

    # Solve variational problem
    solve[a == L, u, bc]

    # Update previous solution
    u_n.assign[u]

5, có nghĩa là tệp

u = Function[V]
t = 0
for n in range[num_steps]:

    # Update current time
    t += dt
    u_D.t = t

    # Solve variational problem
    solve[a == L, u, bc]

    # Update previous solution
    u_n.assign[u]

4 sẽ là một tệp nhỏ duy nhất trỏ đến một số lượng lớn tệp

u = Function[V]
t = 0
for n in range[num_steps]:

    # Update current time
    t += dt
    u_D.t = t

    # Solve variational problem
    solve[a == L, u, bc]

    # Update previous solution
    u_n.assign[u]

5 chứa dữ liệu thực tế. Lưu ý rằng chúng tôi chọn lưu trữ giải pháp vào thư mục con có tên

u_e = interpolate[u_D, V]
error = np.abs[u_e.vector[].array[] - u.vector[].array[]].max[]
print['t = %.2f: error = %.3g' % [t, error]]

0. Điều này là để tránh làm lộn xộn thư mục nguồn của chúng tôi với tất cả các tệp dữ liệu được tạo. Người ta không cần tạo thư mục trước khi chạy chương trình vì nó sẽ được tạo tự động bởi FEniCS

Chương trình hoàn chỉnh xuất hiện bên dưới

from fenics import *
import time

T = 2.0            # final time
num_steps = 50     # number of time steps
dt = T / num_steps # time step size

# Create mesh and define function space
nx = ny = 30
mesh = RectangleMesh[Point[-2, -2], Point[2, 2], nx, ny]
V = FunctionSpace[mesh, 'P', 1]

# Define boundary condition
def boundary[x, on_boundary]:
    return on_boundary

bc = DirichletBC[V, Constant[0], boundary]

# Define initial value
u_0 = Expression['exp[-a*pow[x[0], 2] - a*pow[x[1], 2]]',
                 degree=2, a=5]
u_n = interpolate[u_0, V]

# Define variational problem
u = TrialFunction[V]
v = TestFunction[V]
f = Constant[0]

F = u*v*dx + dt*dot[grad[u], grad[v]]*dx - [u_n + dt*f]*v*dx
a, L = lhs[F], rhs[F]

# Create VTK file for saving solution
vtkfile = File['heat_gaussian/solution.pvd']

# Time-stepping
u = Function[V]
t = 0
for n in range[num_steps]:

    # Update current time
    t += dt

    # Compute solution
    solve[a == L, u, bc]

    # Save to file and plot solution
    vtkfile  0, u[x,t] = 0, x ∈ ∂D, [7] u[x,0] = f [x . Vấn đề sóng 3D là. utt = ∇ 2 u, x ∈ D, t > 0, u[x,t] = 0, x ∈ ∂D, [8] u[x,0] = f [x], x ∈ D, ut [x . . The 3D wave problem is. utt = ∇ 2 u, x ∈ D, t > 0, u[x,t] = 0, x ∈ ∂D, [8] u[x,0] = f [x], x ∈ D, ut [x,0] = g [x], x ∈ D.