# Phương pháp số trong mô phỏng thủy động lực

Thông tin về môn học:
* Bậc học: Cao học
* Thời lượng: 2 tín chỉ
* Thời khóa biểu: từ 20/4/2020
* Hình thức học: online, chia thành 10 bài
Học viên tự bố trí thời gian học phù hợp, thực hiện bài tập.
* Hình thức trao đổi:
* qua Zalo (thông báo chung),
* qua HackMD (trang web này - thông tin về bài học)
* Hình thức đánh giá: bài thu hoạch
* Giảng viên: Nguyễn Quang Chiến, trường ĐH Thủy lợi
* Giáo trình: http://coastalstudy.info/num.meth/book.html
* Tài liệu: http://coastalstudy.info/num.meth/notebooks
(được cập nhật thường xuyên)
* Công cụ học tập: lập trình Julia trực tuyến trên nền CoCalc
https://cocalc.com/

## Bài 1

### Khái niệm chung liên quan đến phương pháp số

Hãy tìm đọc các nội dung sau trong [quyển bài giảng](http://coastalstudy.info/num.meth/book.html).
* Khái niệm về phương pháp số
* Phân biệt phương pháp số và mô hình toán
* Các phương pháp xấp xỉ trong phương pháp số
* Các loại sai số

###### Câu hỏi liên quan đến bài thu hoạch
> ✍ Theo học viên, việc nghiên cứu phương pháp số sẽ giúp ích gì cho việc ứng dụng các mô hình số trị mà học viên từng áp dụng?
>

### Thực hành (1.1)

* Vào trang thực hành trực tuyến [CoCalc](https://cocalc.com/), đăng nhập bằng Google hoặc một cách khác
* Tạo một Project có tên là `Phuong_phap_so` chẳng hạn
* Đưa các file bài mẫu từ [trang web tài liệu](http://coastalstudy.info/num.meth) vào Cocalc. Làm thế nào?
* Trong CoCalc đã có sẵn tính năng Upload file. File có định dạng phù hợp với đuôi `*.ipynb`
* Vì vậy cần tải file `*.ipynb` từ trang web tài liệu về máy của bạn rồi mới upload.
* Thử nghiệm mở một file `*.ipynb` rồi quan sát giao diện trên trang web:
* Tài liệu được mở ra có nhiều ô (cell). Mỗi ô có thể là mã lệnh, kết quả tính, hoặc lời thuyết minh.
* Các ô có thể được Cut - Copy - Paste, thực hiện (Shift+Enter).
* Môi trường CoCalc có thể dùng để thử nghiệm, phù hợp với các chương trình nhỏ, tính toán một chiều.

### Thực hành (1.2)

* Vào trang web http://docs.junolab.org/latest/man/installation/ để cài đặt Julia và môi trường phát triển mã lệnh Juno.
* Môi trường Juno được cài trên máy tính cá nhân của học viên cho phép chạy các chương trình lớn hơn.

###### Ví dụ: Tính lặp chiều dài sóng
> * Làm [bài tập](http://coastalstudy.info/num.meth/book.html#ph%C6%B0%C6%A1ng-tr%C3%ACnh-vi-ph%C3%A2n-th%C6%B0%E1%BB%9Dng).
> * Xem giải thích kĩ hơn và lời giải tại [đường link](https://nbviewer.jupyter.org/url/162.250.127.146/num.meth/notebooks/Dispersion_L.ipynb).
> * Ấn nút tải về file (ở góc trang web) và ấn Ctrl+S để lưu vào máy. Sau đó Upload lên project trong [CoCalc](https://cocalc.com/).
> * Nếu cần, chọn `Julia 1.4` để tính toán theo gợi ý của CoCalc. Đây là phiên bản Julia mới nhất được dùng online.
> * Trong mỗi bảng tính có nhiều ô (Cell). Có hai loại ô:
> * ô mã lệnh (code) viết bằng ngôn ngữ Julia
> * ô chú thích (viết bằng định dạng Markdown)
> Để sửa ô nào ta kích đúp chuột vào ô đó. Cuối cùng, gõ Shift+Enter để hoàn tất. Nếu đó là ô lệnh, câu lệnh sẽ được thực thi.
> * Thực hành thêm vào, xóa, sao chép, di chuyển các ô.

> **Nhàm chán giờ học phương pháp số?**
> Chuyện xảy ra vào tầm năm 2008 khi ĐH Thủy lợi mới thử nghiệm forum của trường để các bạn tâm sự phản ánh cuộc sống kí túc và giảng được. Khi dò qua các bài viết, GV thấy một *thread* có chủ đề như thế này. Hi vọng là đến bây giờ, sau 12 năm, với nền tảng học trực tuyến tốt hơn, với các bài học có ví dụ tính toán cụ thể, sẽ làm cho môn học đỡ lý thuyết nặng nề, "nhàm chán" :D

## Bài 2: Phương trình vi phân thường

Đọc mục [Vi phân số trị](http://coastalstudy.info/num.meth/book.html#div-idsubnum-intvi-ph%C3%A2n-s%E1%BB%91-tr%E1%BB%8B) và làm bài tập sau (viết ra giấy).
> **Bài tập.** ✍ Bằng cách kết hợp các phương trình khai triển Taylor, hãy thiết lập các công thức sai phân lùi và sai phân trung tâm. Công thức sai phân trung tâm có độ chính xác cấp mấy?
> * Sai phân tiến: $$ f'(x) = \frac{f(x+\Delta x) - f(x) } {\Delta x} + O(\Delta x) $$
> * Sai phân lùi: $$ f'(x) = ? $$
> * Sai phân trung tâm: $$ f'(x) = ? $$
>
Cần thấy rằng phương pháp sai phân trung tâm thường có độ chính xác cao hơn sai phân tiến/lùi. Dĩ nhiên, vẫn có những ngoại lệ, song với cùng một hàm $f$ và cùng một độ dài bước $\Delta x$, sai phân trung tâm cho phép ta xấp xỉ đạo hàm chính xác hơn.

> Chú ý:
> * công thức giải tích với hàm liên tục → vi phân
> * công thức xấp xỉ với hàm số tại các điểm rời rạc → sai phân.

###### Bài tập: So sánh các cách tính sai phân
> Hãy tải về [file này](https://nbviewer.jupyter.org/url/162.250.127.146/num.meth/notebooks/Sai_phan.ipynb), tải lên CoCalc rồi thực hành.

> **Ghi chú:** Đối với các phép sai phân có nhiều điểm thì thường được dùng kí hiệu $f_j, f_{j-1}, f_{j+1}$ v.v. Nếu coi lưới điểm có khoảng cách đều thì
> $f_j = f(x), f_{j-1} = f(x - \Delta x), f_{j+1} = f(x + \Delta x).$ Đạo hàm được viết tại điểm $x$, kí hiệu là $f'_j$. Biểu thức tính $f'_j$ còn được gọi là một *sơ đồ sai phân*. Như vậy có sơ đồ sai phân tiến, sơ đồ sai phân lùi, sơ đồ sai phân trung tâm.

* Sơ đồ sai phân tiến: $f_j = (f_{j+1} - f_j) / (\Delta x)$
* Sơ đồ sai phần lùi: $f_j = ?$
* Sơ đồ sai phân trung tâm: $f_j = (f_{j+1} - f_{j-1}) / (2 \Delta x)$

Một khi đã quen với các cách sai phân tiến - lùi - trung tâm, làm bài tập tiếp sau đó:
> **Bài tập.** Rút ra công thức sai phân loại "(1-8-8-1/12)".
>

$$ f'_j = (f_{j+1} - f_{j-1}) / (2 \Delta x) $$

$$ f'_{j+1} = (f_{j+2} - f_{j}) / (2 \Delta x) $$

$$ f'_{j-1} = (f_{j} - f_{j-2}) / (2 \Delta x) $$

Đạo hàm bậc cao: để viết được các biểu thức sai phân bậc hai, cần biến đổi phương trình Taylor.

###### **Bài tập.** Đạo hàm bậc hai
> ✍ Làm thế nào viết được biểu thức tính $f''(x)$ dựa vào hai phương trình Taylor [này](http://coastalstudy.info/num.meth/book.html#eq:taylor1)?
>
> Đó là cách "chính thống". Ngoài ra, còn một cách "mẹo mực" khác, biến đổi gián tiếp qua đạo hàm bậc nhất. Hãy xem quyển [giáo trình](http://coastalstudy.info/num.meth/book.html#th%E1%BB%91ng-nh%E1%BA%A5t-%E1%BB%95n-%C4%91%E1%BB%8Bnh-v%C3%A0-h%E1%BB%99i-t%E1%BB%A5)
>
> 💻 Hãy bổ sung vào trong file Sai_phan một hàm tên là SP2 để tính toán đạo hàm bậc 2, sau đó thử với hai trường hợp có sẵn. Kiểm tra kết quả với:
> * Trường hợp 1: $f''(0) = 2$
> * Trường hợp 2: $f''(0) = 1$

### Khái niệm
* Nắm được khái niệm "sai số cắt cụt" (TE). ✍ Viết sai số cắt cụt cho các công thức sai phân.
* Các khái niệm hội tụ - thống nhất - ổn định: nói chung chỉ để biết rằng mình viết được sơ đồ **hội tụ** là được.

Ví dụ [@Ramshaw2011](http://coastalstudy.info/num.meth/book.html#th%E1%BB%91ng-nh%E1%BA%A5t-%E1%BB%95n-%C4%91%E1%BB%8Bnh-v%C3%A0-h%E1%BB%99i-t%E1%BB%A5) minh họa một trường hợp khi xây dựng công thức sai phân theo cảm tính mà không dựa vào chuỗi Taylor, sẽ dẫn đến công thức sai phân không hội tụ. Ngay lúc này, bạn có thể thấy ví dụ này khó hiểu. Nếu vậy, hãy để ví dụ này lại về sau và đọc tiếp.

### Ngoài sai phân hữu hạn?
* Trong khuôn khổ khóa học này chỉ còn phương pháp *Thể tích hữu hạn* (Thuật ngữ *finite volume* cũng có thể gọi là "khối hữu hạn" mà không khác biệt về ý nghĩa.). Hình vẽ ý tưởng *Thể tích hữu hạn* thể hiện ở [hình vẽ](http://coastalstudy.info/bookfig/pps/FVM_Durran.png).
* Chi tiết xem Mục [Thể tích hữu hạn](http://coastalstudy.info/num.meth/book.html#div-idsecfvm-ph%C6%B0%C6%A1ng-ph%C3%A1p-th%E1%BB%83-t%C3%ADch-h%E1%BB%AFu-h%E1%BA%A1n).

## Bài 3: Hai loại bài toán PT vi phân thường

### Bài toán điều kiện đầu (bài toán dự báo)


Hình trên là minh họa cho bài toán điều kiện đầu. Ta biết tung độ điểm đầu tiên (A~0~). Cần vạch ra đường xanh lam (hàm số). Nhưng ta biết đạo hàm (tức là độ dốc) ứng với giá trị bất kì của *t*.

Dưới dạng toán học: giải PT vi phân

$$ \frac{dy}{dt} = f(t) $$

$$

Biết dạng của hàm f, hãy tìm hàm $y$. Cho trước điều kiện đầu y~0~ tại t~0~.

(Để gắn với ý nghĩa "dự báo" của bài toán này, mà biến độc lập ở đây là *t* thay vì *x*.)

Lấy sai phân PT trên có thể viết được:

$$\frac{y^{n+1} - y^n}{\Delta t} = f(t^n)$$

hay cách khác:

$$ \frac{y^{n+1} - y^n}{\Delta t} = f(t^{n+1}) $$

$$
tùy theo vế phải được chọn theo điểm *t* nào.

#### Phương pháp giải

Phương pháp đơn giản nhất là phương pháp Euler (tiến).

$$ y^{n+1} = y^n + \Delta t f(t^{n})$$

$$

Trên hình minh họa phương pháp Euler tiến. Ở bước thứ nhất, đoạn A~0~ A~1~ được xác định theo hướng tiếp tuyến tại A~0~ và tiếp theo, đoạn A~1~ A~2~ được xác định theo hướng tiếp tuyến tại A~1~ (mũi tên xanh lục).

Nghiệm số trị sẽ là đường A~0~ A~1~ A~2~...

Hiển nhiên, khi co ngắn bước tính (Δt) thì độ chính xác tăng lên. Ngược lại khi Δt ↗ thì càng kém chính xác. Thậm chí, khi Δt tăng đến mức độ nhất định thì nghiệm số trị trở nên bất ổn định.

* Tiêu chí để đánh giá ổn định là hệ số khuếch đại G
* Điều kiện ổn định là |G| < 1
* Biểu hiện của sự bất ổn định nghiệm là đường gấp khúc dao động mạnh và không có xu hướng bắt theo đường cong trơn của nghiệm giải tích.

> &#x1F4BB; **Phương pháp Euler.** Tìm hiểu phương pháp tính Euler từ file [Euler](https://nbviewer.jupyter.org/url/162.250.127.146/num.meth/notebooks/Euler.ipynb). Sau đó áp dụng để giải PT vi phân với trường hợp vế phải là một đa thức trong file này.

Nhận thấy phương pháp Euler có nhược điểm với các hàm mà đồ thị cong; như ở đây ta thấy ngay cả tọa độ A~1~ đã kém chính xác vì được phóng từ A~0~ với độ dốc lớn quá. Nếu độ dốc được chỉnh với thông tin độ dốc tại A~1~ (mũi tên xanh) thì tốt hơn.


Từ đó đưa ra các phương pháp cải thiện khác

> &#x1F4BB; **Cải thiện phương pháp Euler tiến.** Xem file [Bai_toan_DK_dau](https://nbviewer.jupyter.org/url/162.250.127.146/num.meth/notebooks/Bai_toan_DK_dau.ipynb). Từ đó thấy dược ưu điểm của phương pháp cải thiện so với phương pháp Euler gốc.

### Bài toán điều kiện biên

So sánh với bài toán điều kiện đầu: xem hình, trong đó:
1. Bài toán ĐK đầu (hình a): động thái của nghiệm phụ thuộc nhiều vào điều kiện đầu. Có ý nghĩa dự báo. Có trục *t* biểu thị một hướng tiến (thời gian), chính là hướng tiến hành giải dần ra nghiệm PT. Có thể giải mãi (dự đoán mãi về tương lai).

2. Bài toán ĐK biên (hình b): hai đầu bị chặn, động thái của nghiệm bị ràng buộc bởi bản chất PT vi phân. Trục *x* bị giới hạn. Có thể dùng mô tả một hệ thống ỏn định, không biến đổi theo thời gian (phân bố ứng suất, phân bố nhiệt trên vật liệu, dòng chảy ổn định trên một đoạn sông).


Hãy đọc kĩ [dạng bài](http://coastalstudy.info/num.meth/book.html#div-idsecbvpdivb%C3%A0i-to%C3%A1n-%C4%91i%E1%BB%81u-ki%E1%BB%87n-bi%C3%AAn-trong-pvt) trong quyển giáo trình. Hình dung như:

* PT mô tả truyền nhiệt trên thanh
* Đầu bên trái nhiệt độ được giữ cố định &rightarrow; đầu bên trái hình (b) là cố định
* Đầu bên phải nhiệt độ và gradient nhiệt có mối liên hệ &rightarrow; đầu bên phải hình (b) có thể dịch lên xuống nhưng "độ dốc" là không đổi.
* Nhiệt độ mọi điểm có liên hệ với nhau. Khác với bài toán ĐK đầu, ở đó trị số tại mỗi điểm chỉ phụ thuộc vào điểm "bên trái" (tức là quá khứ).

Do sự liên hệ mật thiết giữa các trị số u mà phải giải hệ phương trình. Ở đây là hệ 3 đường chéo. Có thuật toán (Thomas) để giải.

PT sai phân:

$$\frac {-u_{j-1} + u_j - u_{j+1} } {(\Delta x)^2} + q u_j - f_j = 0$$

tương đương với

$$u_{j - 1} -(2 + q (\Delta x)^2) u_j + u_{j+1} = -f_j (\Delta x)^2$$

Đặt $p = -(2 + q (\Delta x)^2)$ và $r_j = -f_j (\Delta x)^2$; hệ PT ba đường chéo như sau: $$\begin{bmatrix}
1 & 0 & & & & \\
1 & p & 1 & & & \\
& 1 & p & 1 & & \\
& & . & . & . & \\
& & & 1 & p & 1 \\
& & & & -1 & 1+\beta \\
\end{bmatrix}
\begin{Bmatrix}
u_1 \\ u_2 \\ u_3 \\ \vdots \\ u_{N-1} \\ u_N
\end{Bmatrix}
=
\begin{Bmatrix}
\alpha \\ r_2 \\ r_3 \\ \vdots \\ r_{N-1} \\ \gamma \Delta x
\end{Bmatrix}$$
$$
> &#x1F4BB; **Phương pháp Thomas.** Thực hành file [Bai_toan_DK_bien](https://nbviewer.jupyter.org/url/162.250.127.146/num.meth/notebooks/Bai_toan_DK_bien.ipynb).

## Bài 4: Phương trình vi phân riêng - PT Elip

Phương trình vi phân riêng (PDE) có chứa các đạo hàm riêng

Lưu ý:

- Bậc của PDE: = cấp đạo hàm cao nhất
- Tuyến tính hoặc phi tuyến:
- Tuyến tính khi hàm *u* và các đạo hàm của nó không nhân với nhau!
- PDE tuyến tính bậc hai có dạng: xem [bài giảng](http://coastalstudy.info/num.meth/book.html#div-idsubclassification-pde-ph%C3%A2n-lo%E1%BA%A1i-pvr)

### Phân loại phương trình vi phân

* Tuyến tính và phi tuyến
* Trong loại tuyến tính: Loại elip, hyperbol và parabol (phân loại theo &Delta;)

Xem kĩ miền ảnh hưởng của từng loại PT vi phân.

### Phương trình "truyền nhiệt" - elip trong không gian 2 chiều

Xem bài giảng ở [Mục này](http://coastalstudy.info/num.meth/book.html#div-idsubdifference-2d-sai-ph%C3%A2n-theo-kh%C3%B4ng-gian-2-chi%E1%BB%81u).

Khi K~x~ và K~y~ không phụ thuộc vào hàm &varphi; thì vấn đề còn lại chỉ còn là sai phân các đạo hàm bậc hai của &varphi; theo *x* và *y*:

$$\frac{\partial \phi}{\partial x^2} = \frac{\phi_{i-1,j} - 2\phi_{i,j} + \phi_{i+1,j}}{(\Delta x)^2} + O(\Delta x^2)$$

$$

###### Bài tập truyền nhiệt trên tấm
> Ví dụ bài toán của [@Chapra2010]. Giải bài toán truyền nhiệt trên tấm vuông cho trước các điều kiện biên như sau:
> 

## Bài 5: Bài toán "chuyển tải" - PT hyperbol

## Bài 6: Sơ đồ sai phân trong không gian 2 chiều

## Bài 7: Phương pháp thể tích hữu hạn

### Thể tích hữu hạn

## Bài 8: Hệ phương trình Navier-Stokes cho chất lỏng không nén

## Bài 9: Hệ phương trình nước nông

## Bài 10: Ôn tập bổ sung phương pháp số cho mô hình 2 chiều