Hướng dẫn newton method optimization python - python tối ưu hóa phương pháp newton

Nguồn: tiểu sử.com

Lời giải thích ngắn nhất và rõ ràng nhất từ ​​đường thẳng

Phương pháp Newton Newton để tối ưu hóa là một trường hợp cụ thể của phương pháp gốc.

Với F F ′ (xk) là đạo hàm của đạo hàm của F F được đánh giá tại Lặp lại.f′′(xk​)” being the derivative of the derivative of “f” evaluated at iteration “k”.

Xem xét chức năng:

Hãy đưa ra âm mưu:

import matplotlib.pyplot as plt
import numpy as np
from numpy import log
x = np.linspace(0, 1.2, 100)
y = 2*x - log(x)
plt.plot(x, y, '-r', label='y = 2*x - log(x)')
plt.title('Graph of y = 2*x - log(x)')
plt.xlabel('x', color='#1C2833')
plt.ylabel('y', color='#1C2833')
plt.legend(loc='upper left')
plt.grid()
plt.show()

Các dẫn xuất thứ nhất và thứ hai là:

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

Hiển thị trong chế độ xem bình thường:

Hãy bắt đầu lặp lại từ x0 = 0,1:

x0 = 0.1fx0 = 2*x0 - ln(x0)
print(fx0)
# 2.50258509299405

Để tìm X1, hãy sử dụng công thức từ trên cùng:

x1 = x0 - (2 - 1/x0) / x0**(-2)
print(x1)
# 0.18000000000000002fx1 = 2*x1 - ln(x1)
print(fx1)
# 2.07479842809193

Để tìm x2 chúng ta nên reepeat nó:

x2 = x1 - (2 - 1/x1) / x1**(-2)
print(x2)
# 0.2952fx2 = 2*x2 - ln(x2)
print(fx2)
# 1.81050218625582

Trông nhàm chán. Cho phép tạo một vòng lặp và kết quả in:

Cho phép kiểm tra chức năng thực tối thiểu:

from sympy.calculus.util import Interval, minimum
from sympy.abc import x
from sympy import ln
print(minimum(2*x - ln(x), x, Interval(0, 1)))
# log(2) + 1

log (2) + 1 = 1.6931471805599454, vì vậy độ chính xác trên lần lặp thứ năm là khoảng 0,0000003.

Hãy thể hiện vẻ ngoài của nó:

for x, fx in zip(x_list, fx_list):
    plt.plot(x, fx, marker="o", markersize=10, markeredgecolor="red", markerfacecolor="green")plt.show()

Mã đầy đủ ở đây:

Nếu bạn cần một số công cụ tối ưu hóa 3D, hãy thử xem:

Phương pháp Newton-Raphson

(Nguồn: https://phvu.net)

Trong chương trình phổ thông, khoảng lớp 5 học sinh được học cách giải phương trình bậc nhất 1 ẩn (qua một lớp bài toán được biết với cái tên nổi tiếng là tìm x). Khoảng cấp 2, HS được học cách giải phương trình bậc hai một ẩn và hệ phương trình bậc nhất hai ẩn.

Tốt nghiệp cấp 3, nếu đi học tiếp Cử nhân Công nghệ thông tin (như người viết bài này), HS sẽ được học tiếp cách giải hệ phương trình tuyến tính nhiều ẩn dạng (với là ma trận ) trong môn Đại số tuyến tính. Giải thuật lúc đó được dạy thì người viết không còn nhớ rõ, nhưng đại khái là tìm cách sử dụng các phép biến đổi dòng trên ma trận một cách khéo léo để tìm nghiệm. Mặc dù phương pháp đó có thể cài đặt được thành thuật toán, nhưng vì dựa vào trick (và “quan sát” theo kiểu của con người) khá nhiều nên không thực sự straightforward (sau vài ngày thì người viết nhớ ra tên phương pháp là Phép khử Gaussian) (với là ma trận ) trong môn Đại số tuyến tính. Giải thuật lúc đó được dạy thì người viết không còn nhớ rõ, nhưng đại khái là tìm cách sử dụng các phép biến đổi dòng trên ma trận một cách khéo léo để tìm nghiệm. Mặc dù phương pháp đó có thể cài đặt được thành thuật toán, nhưng vì dựa vào trick (và “quan sát” theo kiểu của con người) khá nhiều nên không thực sự straightforward (sau vài ngày thì người viết nhớ ra tên phương pháp là Phép khử Gaussian)

Người viết có dịp tìm hiểu phương pháp Newton để giải hệ phương trình phi tuyến lần đầu tiên trong quá trình thực hiện luận văn. Sau đó gần đây lại có cơ hội làm việc với các phương pháp khác như Conjugate Gradient, LSQR …

Một điều đáng ngạc nhiên là, theo quá trình học của bản thân người viết, những phương pháp giải xấp xỉ/chính xác các phương trình/hệ phương trình tuyến tính/phi tuyến đã không được hệ thống hoá một cách đầy đủ, cả trong chương trình Toán phổ thông lẫn Đại học không chuyên ngành Toán (mặc dù tư tưởng chung của các phương pháp này khá tự nhiên). Vì vậy chuỗi bài này muốn hệ thống lại các phương pháp trên. Bài đầu tiên sẽ nói về phương pháp Newton-Raphson. Các thuật toán khác và một số vấn đề liên quan sẽ được trình bày trong các bài sau. Các phương pháp này là bước đi ban đầu vào optimization, lĩnh vực mà người viết mong muốn được tìm hiểu sâu hơn.

Phương pháp Newton-Raphson giải một (hệ) phương trình thực bằng cách xấp xỉ dần giá trị của nghiệm thông qua một dãy số (hội tụ về giá trị thật của nghiệm?). Chúng ta sẽ bắt đầu với trường hợp hàm một biến, sau đó mở rộng cho trường hợp nhiều biến.

1. Trường hợp một biến

Cho hàm một biến , mục tiêu là tìm giá trị của thỏa phương trình . (Tất nhiên trong trường hợp là hàm bậc nhất theo x thì không còn gì phải bàn, ở đây ta muốn xét đến trường hợp tổng quát)., mục tiêu là tìm giá trị của

Phương pháp Newton xuất phát từ việc xấp xỉ giá trị của đạo hàm tại một điểm:

Do ta đang muốn tìm nghiệm sao cho , do đó với thì phương trình trên trở thành: sao cho , do đó với thì phương trình trên trở thành:

Tổng quát, phương pháp Newton lặp liên tục để tìm nghiệm ngày càng gần giá trị thực, với công thức nghiệm là:

với giá trị khởi tạo nào đó. nào đó.

Hình 1 cho một cái nhìn trực quan hơn về lời giải của phương pháp này. Với càng lớn, giá trị của càng gần về vị trí giao với trục hoành của đồ thị hàm . càng lớn, giá trị của càng gần về vị trí giao với trục hoành của đồ thị hàm .

Tuy nhiên cũng giống như nhiều phương pháp tính xấp xỉ khác, khả năng tốc độ hội tụ của phương pháp Newton phụ thuộc nhiều vào giá trị “đoán mò” ban đầu . Giá trị càng gần nghiệm thực của phương trình thì phương pháp này càng hội tụ nhanh chóng. Trong trường hợp xấu (khi vô nghiệm, hoặc khi giá trị đoán mò không cùng nằm trên một đoạn biến thiên của hàm số) thì phương pháp sẽ khó hội tụ nhanh. Trong trường hợp phương trình không có nghiệm, phương pháp này cũng không thể phát hiện.. Giá trị càng gần nghiệm thực của phương trình thì phương pháp này càng hội tụ nhanh chóng. Trong trường hợp xấu (khi vô nghiệm, hoặc khi giá trị đoán mò không cùng nằm trên một đoạn biến thiên của hàm số) thì phương pháp sẽ khó hội tụ nhanh. Trong trường hợp phương trình không có nghiệm, phương pháp này cũng không thể phát hiện.

Một điểm hạn chế nữa là do yêu cầu tính tích phân của hàm số, phương pháp này sẽ khó khăn trong trường hợp việc tính đạo hàm trở nên phức tạp. Trong trường hợp này, người ta có thể xấp xỉ việc tính đạo hàm bằng cách tìm phương trình đường thẳng đi qua hai điểm lân cận của hàm số tại vị trí đang xét (để xấp xỉ hệ số góc của tiếp tuyến – vốn có thể tính chính xác bằng đạo hàm tại điểm đó).

Do sự phức tạp trong việc tính đạo hàm nên việc cài đặt một chương trình tổng quát cho thuật toán Newton-Raphson là tương đối khó. Tuy nhiên trong trường hợp là đa thức theo thì có thể dễ dàng cài đặt phương pháp này trong MATLAB như sau là đa thức theo thì có thể dễ dàng cài đặt phương pháp này trong MATLAB như sau

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

x0 = 0.1fx0 = 2*x0 - ln(x0)
print(fx0)
# 2.50258509299405

x0 = 0.1fx0 = 2*x0 - ln(x0)
print(fx0)
# 2.50258509299405

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

x0 = 0.1fx0 = 2*x0 - ln(x0)
print(fx0)
# 2.50258509299405

x0 = 0.1fx0 = 2*x0 - ln(x0)
print(fx0)
# 2.50258509299405

x0 = 0.1fx0 = 2*x0 - ln(x0)
print(fx0)
# 2.50258509299405

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

x0 = 0.1fx0 = 2*x0 - ln(x0)
print(fx0)
# 2.50258509299405

x0 = 0.1fx0 = 2*x0 - ln(x0)
print(fx0)
# 2.50258509299405

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

x1 = x0 - (2 - 1/x0) / x0**(-2)
print(x1)
# 0.18000000000000002fx1 = 2*x1 - ln(x1)
print(fx1)
# 2.07479842809193

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

x1 = x0 - (2 - 1/x0) / x0**(-2)
print(x1)
# 0.18000000000000002fx1 = 2*x1 - ln(x1)
print(fx1)
# 2.07479842809193

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

x1 = x0 - (2 - 1/x0) / x0**(-2)
print(x1)
# 0.18000000000000002fx1 = 2*x1 - ln(x1)
print(fx1)
# 2.07479842809193

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

x1 = x0 - (2 - 1/x0) / x0**(-2)
print(x1)
# 0.18000000000000002fx1 = 2*x1 - ln(x1)
print(fx1)
# 2.07479842809193

x1 = x0 - (2 - 1/x0) / x0**(-2)
print(x1)
# 0.18000000000000002fx1 = 2*x1 - ln(x1)
print(fx1)
# 2.07479842809193

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

x2 = x1 - (2 - 1/x1) / x1**(-2)
print(x2)
# 0.2952fx2 = 2*x2 - ln(x2)
print(fx2)
# 1.81050218625582

x2 = x1 - (2 - 1/x1) / x1**(-2)
print(x2)
# 0.2952fx2 = 2*x2 - ln(x2)
print(fx2)
# 1.81050218625582

x2 = x1 - (2 - 1/x1) / x1**(-2)
print(x2)
# 0.2952fx2 = 2*x2 - ln(x2)
print(fx2)
# 1.81050218625582

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

x1 = x0 - (2 - 1/x0) / x0**(-2)
print(x1)
# 0.18000000000000002fx1 = 2*x1 - ln(x1)
print(fx1)
# 2.07479842809193

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

x1 = x0 - (2 - 1/x0) / x0**(-2)
print(x1)
# 0.18000000000000002fx1 = 2*x1 - ln(x1)
print(fx1)
# 2.07479842809193

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

x2 = x1 - (2 - 1/x1) / x1**(-2)
print(x2)
# 0.2952fx2 = 2*x2 - ln(x2)
print(fx2)
# 1.81050218625582

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

from sympy.calculus.util import Interval, minimum
from sympy.abc import x
from sympy import ln
print(minimum(2*x - ln(x), x, Interval(0, 1)))
# log(2) + 1

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

from sympy.calculus.util import Interval, minimum
from sympy.abc import x
from sympy import ln
print(minimum(2*x - ln(x), x, Interval(0, 1)))
# log(2) + 1

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

from sympy.calculus.util import Interval, minimum
from sympy.abc import x
from sympy import ln
print(minimum(2*x - ln(x), x, Interval(0, 1)))
# log(2) + 1

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

from sympy.calculus.util import Interval, minimum
from sympy.abc import x
from sympy import ln
print(minimum(2*x - ln(x), x, Interval(0, 1)))
# log(2) + 1

from sympy.calculus.util import Interval, minimum
from sympy.abc import x
from sympy import ln
print(minimum(2*x - ln(x), x, Interval(0, 1)))
# log(2) + 1

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

from sympy.calculus.util import Interval, minimum
from sympy.abc import x
from sympy import ln
print(minimum(2*x - ln(x), x, Interval(0, 1)))
# log(2) + 1

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

from sympy.calculus.util import Interval, minimum
from sympy.abc import x
from sympy import ln
print(minimum(2*x - ln(x), x, Interval(0, 1)))
# log(2) + 1

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

for x, fx in zip(x_list, fx_list):
    plt.plot(x, fx, marker="o", markersize=10, markeredgecolor="red", markerfacecolor="green")plt.show()

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

x1 = x0 - (2 - 1/x0) / x0**(-2)
print(x1)
# 0.18000000000000002fx1 = 2*x1 - ln(x1)
print(fx1)
# 2.07479842809193

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

for x, fx in zip(x_list, fx_list):
    plt.plot(x, fx, marker="o", markersize=10, markeredgecolor="red", markerfacecolor="green")plt.show()

x1 = x0 - (2 - 1/x0) / x0**(-2)
print(x1)
# 0.18000000000000002fx1 = 2*x1 - ln(x1)
print(fx1)
# 2.07479842809193

1. Đoán mò nghiệm 

2. Tính ma trận Jacobi và vector thặng dư.

3. Giải phương trình (2) bằng các phương pháp trong đại số tuyến tính

4. Cập nhật giá trị  theo (6).

5. Lặp lại bước 2 nếu lời giải chưa hội tụ.

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

1. Đoán mò nghiệm 

2. Tính ma trận Jacobi và vector thặng dư.

3. Giải phương trình (2) bằng các phương pháp trong đại số tuyến tính

4. Cập nhật giá trị  theo (6).

5. Lặp lại bước 2 nếu lời giải chưa hội tụ.

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

x1 = x0 - (2 - 1/x0) / x0**(-2)
print(x1)
# 0.18000000000000002fx1 = 2*x1 - ln(x1)
print(fx1)
# 2.07479842809193

1. Đoán mò nghiệm 

2. Tính ma trận Jacobi và vector thặng dư.

3. Giải phương trình (2) bằng các phương pháp trong đại số tuyến tính

4. Cập nhật giá trị  theo (6).

5. Lặp lại bước 2 nếu lời giải chưa hội tụ.

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

1. Đoán mò nghiệm 

2. Tính ma trận Jacobi và vector thặng dư.

3. Giải phương trình (2) bằng các phương pháp trong đại số tuyến tính

4. Cập nhật giá trị  theo (6).

5. Lặp lại bước 2 nếu lời giải chưa hội tụ.

x2 = x1 - (2 - 1/x1) / x1**(-2)
print(x2)
# 0.2952fx2 = 2*x2 - ln(x2)
print(fx2)
# 1.81050218625582

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

x1 = x0 - (2 - 1/x0) / x0**(-2)
print(x1)
# 0.18000000000000002fx1 = 2*x1 - ln(x1)
print(fx1)
# 2.07479842809193

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

x1 = x0 - (2 - 1/x0) / x0**(-2)
print(x1)
# 0.18000000000000002fx1 = 2*x1 - ln(x1)
print(fx1)
# 2.07479842809193

x1 = x0 - (2 - 1/x0) / x0**(-2)
print(x1)
# 0.18000000000000002fx1 = 2*x1 - ln(x1)
print(fx1)
# 2.07479842809193

Trong chương trình trên, là vector chứa hệ số của đa thức. Chẳng hạn với đa thức thì . Tham số là độ lỗi cho phép, chương trình dừng khi sai số nhỏ hơn ngưỡng này. Trong trường hợp không tìm được nghiệm, chương trình dừng sau lần lặp. là vector chứa hệ số của đa thức. Chẳng hạn với đa thức thì . Tham số là độ lỗi cho phép, chương trình dừng khi sai số nhỏ hơn ngưỡng này. Trong trường hợp không tìm được nghiệm, chương trình dừng sau lần lặp.

Sử dụng hàm này, ta có thể tính nhanh nghiệm của phương trình (cũng chính là tìm xấp xỉ giá trị ) như sau: (cũng chính là tìm xấp xỉ giá trị ) như sau:

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

Một điều khá thú vị là sử dụng cặp hàm , ta thấy hàm chạy không chậm hơn phép lũy thừa trong MATLAB nhiều lắm khi tính (xét tới việc ta đã “đoán mò” giá trị rất thô thiển)., ta thấy hàm chạy không chậm hơn phép lũy thừa trong MATLAB nhiều lắm khi tính (xét tới việc ta đã “đoán mò” giá trị rất thô thiển).

from sympy import diff, ln
from sympy.abc import x

f = 2*x - ln(x)
df = diff(f)
ddf = diff(df)
print(df)
print(ddf)
# 2 - 1/x
# x**(-2)

2. Trường hợp nhiều biến

Giả sử ta phải giải hệ phương trình biến sau: phương trình biến sau:

Một trong những phương pháp giải hệ phương trình này là phương pháp Newton nhiều biến (Multivariate Newton-Raphson Method – MNRM). Mở rộng từ phương pháp Newton cho hàm một biến, ta cũng bắt đầu từ việc đạo hàm từng phần mỗi phương trình ::

Giống như trước, ta xấp xỉ như sau:

              (1)(1)

với là giá trị của nghiệm ở lần lặp thứ k. Ta mong muốn , do đó hệ phương trình (1) có thể viết dưới dạng ma trận: là giá trị của nghiệm ở lần lặp thứ k. Ta mong muốn , do đó hệ phương trình (1) có thể viết dưới dạng ma trận:

      (2)(2)

trong đó là ma trận Jacobi kích thước tại bước lặp thứ là ma trận Jacobi kích thước tại bước lặp thứ

      (3)(3)

là vector thặng dư kích thước tại bước lặp thứ : tại bước lặp thứ :

        (4)(4)

và cuối cùng

            (5)(5)

Do đó nghiệm của hệ tại bước thứ là là

         (6)(6)

Sử dụng những phương trình này, thuật toán Newton tìm lời giải cho hệ phương trình rất đơn giản:

1. Đoán mò nghiệm 

2. Tính ma trận Jacobi và vector thặng dư.

3. Giải phương trình (2) bằng các phương pháp trong đại số tuyến tính

4. Cập nhật giá trị  theo (6).

5. Lặp lại bước 2 nếu lời giải chưa hội tụ.

3. Tổng kết

Phương pháp Newton-Raphson trong bài là một phương pháp hết sức đơn giản. Nó nằm trong họ các phương pháp Household mà nếu có dịp sẽ quay trở lại sau. Do hạn chế trong việc phải tính được đạo hàm của hàm số, phương pháp này thực sự có ý nghĩa khi biết trước dạng hàm. Trong trường hợp không biết trước dạng hàm thì cũng có thể áp dụng phương pháp này bằng cách tính xấp xỉ giá trị đạo hàm tại từng điểm (một trong những cách xấp xỉ là tính hiệu với “đủ gần” với ). Phương pháp xấp xỉ đạo hàm này đặc biệt có ý nghĩa khi cài đặt trên máy tính. với “đủ gần” với ). Phương pháp xấp xỉ đạo hàm này đặc biệt có ý nghĩa khi cài đặt trên máy tính.

Trong MATLAB, việc tìm nghiệm của hệ các phương trình được cài đặt sẵn trong hàm fslove(), fzero()… Các hàm này cài đặt một số thuật toán khác để không phải tính đạo hàm như thuật toán Newton, hoặc tìm cách xấp xỉ giá trị đạo hàm. Đặc biệt hàm fsolve() không cần giá trị đoán mò ban đầu.

Một trường hợp đơn giản là giải hệ phương trình tuyến tính bậc nhất thì trong MATLAB được cài bằng phép chia trái , một số thư viện số học khác như LAPACK, gsl, BLAS cũng có sẵn hàm cho trường hợp này. Trường hợp có nhiều phương pháp giải, trong đó có Conjugate gradient, Cholesky decomposition… Những phương pháp này sẽ là chủ đề của bài tiếp theo. thì trong MATLAB được cài bằng phép chia trái , một số thư viện số học khác như LAPACK, gsl, BLAS cũng có sẵn hàm cho trường hợp này. Trường hợp có nhiều phương pháp giải, trong đó có Conjugate gradient, Cholesky decomposition… Những phương pháp này sẽ là chủ đề của bài tiếp theo.

———–&&———–

Filed under: Giải tích số |