Bài giảng Chuyên đề Phương pháp tính - Chương 6: Nghiệm gần đúng của hệ phương trình vi phân thường

Khoa Xây Dựng Thủy Lợi Thủy Điện

Bộ môn Cơ Sở Kỹ Thuật

NGHIỆM GẦN ĐÚNG CỦA HỆ

^{Chương 6}PHƯƠNG TRÌNH VI PHÂN THƯỜNG

6.1 Mở đầu

Nhiều bài toán khoa học kỹ thuật có phương trình chỉ đạo là (hệ)

phương trình vi phân thường cùng với điều kiện biên và điều kiện ban đầu.

Nghiệm đúng của chúng thường chỉ áp dụng cho một số lớp bài toán rất hạn

chế; đa số các bài toán là phải tìm nghiệm gần đúng.

Có hai loại bài toán là:

(i)

Bài toán Cauchy hay còn gọi là bài toán giá trị ban đầu, bao gồm

(hệ) phương trình vi phân và điền kiện ban đầu của bài toán.

(ii) Bài toán biên, bao gồm (hệ) phương trình vi phân và điều kiện biên

Để giải gần đúng các bài toán nầy có hai phương pháp là:

(a) Phương pháp giải tích: tìm nghiệm gần đúng dưới dạng biểu thức

như phương pháp xấp xỉ liên tiếp Picard, phương pháp chuổi

nguyên, phương pháp tham số bé,…

(b) Phương pháp số: tìm nghiệm gần đúng bằng số tại các điểm rời rạc;

nó còn chia ra phương pháp một bước (như phương pháp Euler,

Runghe-Kutta,…) và phương pháp đa bước (Adams,…); Với

phương pháp một bước tính nghiệm gần đúng y_ithông qua y_i-1còn

với phương pháp đa bước y_itính được thông qua nhiều bước trước

đó: y_i-1, y_i-2, y_i-3,…

6.2 Nghiệm gần đúng của bài toán Cauchy đối với phương trình vi

phân thường

y'= f(x, y)

y(x₀) = y₀







Giả sử ta cần giải bài toán Cauchy:

(6.2.1)

Giả sử rằng trong miền ta xét, hàm f(x,y) có các đạo hàm riêng liên

tục đến cấp n, khi đó nghiệm cần tìm sẽ có các đạo hàm riêng liên tục đến

cấp n + 1, và do đó ta có thể viết :

∆y = y(x ) − y = (x − x )y^,+

o

0

(x − x₀)²

(x − x₀)ⁿ⁺¹

(n +1)!

n+1

y₀⁽ⁿ⁺¹⁾+θ( x − x₀

)

(6.2.2)

y"₀+...... +

2!

Ký hiệu x - x₀= h, với h đủ bé ta có thể bỏ qua 0(|x – x₀|ⁿ⁺¹).

n+1

θ( x − x₀)

Bài Giảng Chuyên Đề Phương Pháp Tính

Trang 33

Khoa Xây Dựng Thủy Lợi Thủy Điện

Bộ môn Cơ Sở Kỹ Thuật

Từ (6.2.2) ta có: ∆y₀= y(x₀+h) - y₀+ hy’₀+

h²

2!

hⁿ⁺¹

(n +1)!

(n+1)

y"₀+..........+

y₀

(6.2.3)

Để tính (6.2.3) ta lần lượt tính từ (6.2.1):

∂f₀

∂x

∂f₀

∂y

+ f₀

y’₀= f(x₀,y₀) = f₀, y”₀=

,

^m

∂

∂^mu

n





+ f  u = C^Kf ^K

Nói chung ta có:

∑

h^K

m

∂x

∂y

K=0



∂x^m−K∂y^K





n

y^(K)(x )

Vậy ta tính được: y(x) ≅

∑

0

K!

Trong thực tế cách tính nầy^Kí⁼t⁰dùng vì cồng kềnh; ta sẽ xét các phương

pháp giải khác đơn giản hơn.

6.2.1 Phương pháp xấp xỉ liên tiếp Pica

Một trong những phương pháp giải tích giải gần đúng phương trình vi

phân (6.2.1) là phương pháp xấp xỉ liên tiếp Pica.

Mục đích của phương pháp này là xây dựng nghiệm cần tìm là y= y(x)

x

^xdy = f(t,y)dt

⇒ y(x) − y(x₀) = f(t,y)dt

Từ (6.2.1) ta có:

∫ ∫

∫

x₀

x

y(x) = y₀+ f (t,y)dt

Hay:

(6.2.4)

∫

x₀

∂f

∂y

Giả sử f(x,y) là hàm liên tục theo x,y và

< K.

Để tìm xấp xỉ liên tiếp, trong (6.2.4) thay y bằng y₀, ta có xấp xỉ thứ

nhất:

x

y = y + f(t,y )dt

,

1

0

∫

x₀

x

y₂= y₀+ f(t,y₁)dt

Tương tự có xấp xỉ thứ hai:

Tổng quát, ta có:

∫

x₀

x

y_n= y₀+ f(t, y_n−1)dt

, với n = 1,2,3,…

∫

x₀

Bài Giảng Chuyên Đề Phương Pháp Tính

Trang 34

Khoa Xây Dựng Thủy Lợi Thủy Điện

Bộ môn Cơ Sở Kỹ Thuật

x

y(x) ≈ y_n(x) = y₀+ f(t, y_n−1)dt

Như vậy ta sẽ có:

∫

x₀

lim y_n(x) = y(x)

n→∞

M(KC)ⁿ

Sai số: y_n(x) − y(x) ≤

, trong đó f(x,y) = M

K.n!

b

M







a,





x − x₀

y − y₀

Với:

< a ≤ ∞,

< b ≤ ∞ , thì C = min



Ta có:

∂f

∂y

∂f

(i)

> 0 và f(x,y₀) > 0 thì: y₀< y₁< y₂< . . . < y_n< y(x)

> 0 và f(x,y₀) < 0 thì: y₀> y₁> y₂> . . . > y_n> y(x)

(ii)

∂y

Trong hai trường hợp nầy ta có dãy xấp xỉ 1 phía.

∂f

(iii)

< 0 các xấp xỉ Pica lập thành các xấp xỉ 2 phía.

∂y

6.2.2 Phương pháp Euler

y

y=f(x)

A3

A2

`

A1

Ao

x

O

x

3

xo

x

1

x

2

Trước hết chia đọan [x_o, X] thành n đọan nhỏ:

x_i=x_o+ih, với i = 0,1,2,....,n

(X − x_o)

h =

n

Đi xây dựng công thức, dùng khai triển Taylor hàm y=f(x) tại x_ita có:

Bài Giảng Chuyên Đề Phương Pháp Tính

Trang 35

Khoa Xây Dựng Thủy Lợi Thủy Điện

Bộ môn Cơ Sở Kỹ Thuật

′′

2!

y (c_i)

(x − x_i)²

′

y(x) = y(x_i) + y (x_i).(x - x_i) +

Với: c_i= x_i+ θ(x - x_i), 0 < θ < 1

Thay x = x_i+1= x_i+ h, và y^’(x_i) = f(x_i,y(x_i))

′′

2!

y (c_i)

y(x_i+1) = y(x_i) + h.f(x_i, y(x_i)) + h².

Ta có:

Khi bước chia h khá bé, số hạng cuối ≅ 0, khi thay y(x_i) bằng u_ita được:

u_i+1= u_i+ h_i.f(x_i,u_i)

Biểu thức nầy cho phép tính u_i+1khi biết u_i,với điều kiện ban đầu được cho

là: u_o= η

Đánh giá sai số:

∂f

Định lý: Gỉa sử

≤ L và y^''≤ K, trong đó L, K là những hằng số, khi

∂y

đó phương pháp Euler hội tụ và sai số là e_i= u_i- y(x_i) có đánh giá:

e_i= u_i− y(x_i) ≤ M(e₀+ αh)

K

2

M = e^{L(x −x )},α =

i

0

6.2.3 Phương pháp Runghe - Kutta bậc 4

Xét phương trình vi phân:

u’ = f(x , u)

k = h.f(x ,u )



1

i





k = h.f(x + 0.5h,u + 0.5k )

1

6

2

i

1

(k₁+ 2k₂+ 2k₃+ k₄)

⇒ u_{i +1}= u_i+

k₃= h.f(x_i+ 0.5h,u_i+ 0.5k₂)

k₄= h.f(x_i+ h,u_i+ k₃)





u_i− Y(x_i) = 0(h⁴)

Với sai số:

6.2.4 Phương pháp Adam

Giả sử cần giải phương trình vi phân:

Y’ = f(x , y), với điều kiện ban đầu: y(x₀) = y₀

Cho biến số thay đổi bởi bước h nào đó; xuất phát từ điều kiện ban

đầu Y(x₀) = Y₀bằng phương pháp nào đó (ví dụ: phương pháp Runghe-

Kutta bậc 4), ta tìm được 3 giá trị tiếp theo của hàm cần tìm y(x): Y₁=

Y(x₁) = Y(x₀+h), Y₂= Y(x₀+2h), Y₃= Y(x₀+ 3h) .

Bài Giảng Chuyên Đề Phương Pháp Tính

Trang 36

Khoa Xây Dựng Thủy Lợi Thủy Điện

Bộ môn Cơ Sở Kỹ Thuật

Nhờ các giá trị x₀, x₁, x₂, x₃và Y₀, Y₁, Y₂, Y₃, ta tính được q₀, q₁,

q₂, q₃.

’

Trong đó: q₀= h.Y₀= h.f(x₀, y₀), q₁= h.f(x₁, y₁), q₂= h.f(x₂, y₂),

q₃= h.f(x₃, y₃), sau đó ta lập bảng sai phân hữu hạn của các đại lượng y và q

x

y

q

-------

∆y

∆y_o

∆y₁

∆y₂

∆q

∆q₀

∆q₁

∆q₂

∆²q

∆³q

x_o

y_o

q_o

∆²q₀

∆²q₁

x₁

x₂

x₃

y₁

q₁

q₂

∆³q₀

y₂

y₃

--------

-------

--------

q₃

--------- --------- --------- --------- --------- --------- --------- -------

---

-

--

-

Biết các số ở đường chéo dưới, ta tìm ∆y₃theo công thức Adam như sau:

1

2

5

12

3

8

∆y₃= q₃+ .∆q₂+ .∆²q₁+ .∆³.q₀

Tiếp đó ta có:

Y₄= Y₃+ ∆Y₃→ q₄= h.f(x₄, Y₄)

Sau đó viết đường chéo tiếp theo như sau:

∆q₃= q₄- q₃, ∆²q₂= ∆q₃- ∆q₂, ∆³q₁= ∆².q₂- ∆².q

Đường chéo mới cho phép ta tính ∆Y₄:

∆Y₄= q₄+ 1/2∆q₃+ 5/12∆²q₂+ 3/8∆³q₁

Vì vậy ta có: Y₅= Y₄+ ∆Y₄. . . . .

Bài Giảng Chuyên Đề Phương Pháp Tính

Trang 37