2011-01-11 47 views
8

Điều này sẽ rất đơn giản. Tôi có chức năng f(x) và tôi muốn đánh giá f'(x) cho một số x nhất định trong MATLAB.cách đánh giá đạo hàm của hàm trong MATLAB?

Tất cả các tìm kiếm của tôi đều có tính toán biểu tượng, không phải là những gì tôi cần, tôi cần sự khác biệt về số.

Ví dụ: nếu tôi xác định: fx = inline('x.^2')

Tôi muốn tìm nói f'(3), đó sẽ là 6, tôi không muốn tìm 2x

Trả lời

7

Để có được một sự khác biệt số (chênh lệch đối xứng), bạn tính toán (f(x+dx)-f(x-dx))/(2*dx)

fx = @(x)x.^2; 
fPrimeAt3 = (fx(3.1)-fx(2.9))/0.2; 

Ngoài ra, bạn có thể tạo ra một vector của các giá trị chức năng và áp dụng DIFF, tức là

xValues = 2:0.1:4; 
fValues = fx(xValues); 
df = diff(fValues)./0.1; 

Lưu ý rằng diff mất sự khác biệt về phía trước và giả định rằng dx bằng 1.

Tuy nhiên, trong trường hợp, bạn có thể tốt hơn để xác định fxpolynomial và đánh giá đạo hàm của hàm, thay vì giá trị hàm.

+0

+1 cho câu trả lời gọn gàng :) – posdef

+0

Tôi có nghĩ rằng tôi chọn dx nhỏ hơn, câu trả lời của tôi càng chính xác hơn. Nếu vậy có một hằng số MATLAB cho một số thực rất nhỏ? Một cái gì đó giống như pi cho 3,14 ... hoặc i cho sqrt (-1)? – lms

+2

@codenoob: 'eps' cung cấp cho bạn một số lượng rất nhỏ. Tuy nhiên, đối với hầu hết các mục đích thực tế, 0,0001 là đủ. Ngoài ra, nếu bạn đang làm việc với đa thức và sử dụng 'polyder', bạn không phải lo lắng về kích thước của' dx'. – Jonas

4

bạn đã thử diff (tính toán sự khác biệt và xấp xỉ một dẫn xuất), gradient, hoặc polyder (tính toán đạo hàm của một hàm đa thức)?

Bạn có thể đọc thêm về các chức năng này bằng cách sử dụng help <commandname> trên bảng điều khiển MATLAB hoặc sử dụng trình duyệt chức năng trong menu Trợ giúp.

+0

+1 để được nhanh hơn một chút :) – Jonas

0

Đối với một chức năng được đưa ra trong hình thức phân tích, bạn có thể đánh giá đạo hàm tại một điểm mong muốn với đoạn mã sau:

syms x 
df = diff(x^2); 
df3 = subs(df, 'x', 3); 
fprintf('f''(3)=%f\n', df3); 

Đối với các dẫn xuất số thuần túy sử dụng các giải pháp đã được đưa ra bởi Jonas và posdef.

6

Thiếu hộp công cụ tượng trưng, ​​không có gì ngăn bạn sử dụng Derivest, một công cụ để phân biệt số tự động thích ứng.

derivest(@sin,pi) 
ans = 
      -1 

Ví dụ của bạn rất độc đáo. Trong thực tế, nó thậm chí cung cấp một ước tính của lỗi trong xấp xỉ kết quả.

fx = inline('x.^2'); 
[fp,errest] = derivest(fx,3) 

fp = 
      6 
errest = 
    3.6308e-14 
+0

Điều này thật thú vị, thx. Tôi đã thực hiện nó bằng cách sử dụng phương pháp khác biệt đối xứng mặc dù. – lms

9

Nếu chức năng của bạn được biết đến là hai lần khả vi, sử dụng

f'(x) = (f(x + h) - f(x - h))/2h 

đó là trật tự thứ hai chính xác trong h. Nếu nó chỉ có thể khác biệt một lần, hãy sử dụng

f'(x) = (f(x + h) - f(x))/h  (*) 

là thứ tự đầu tiên trong h.

Đây là lý thuyết. Trong thực tế, mọi thứ khá phức tạp. Tôi sẽ lấy công thức thứ hai (thứ tự đầu tiên) khi phân tích đơn giản hơn. Làm thứ tự thứ hai như một bài tập.

Quan sát đầu tiên là bạn phải đảm bảo rằng (x + h) - x = h, nếu không bạn sẽ gặp phải các lỗi rất lớn. Thật vậy, f (x + h) và f (x) gần nhau (2.0456 và 2.0467), và khi bạn trừ chúng, bạn sẽ mất rất nhiều số liệu quan trọng (ở đây là 0,0011, có 3 số liệu quan trọng ít hơn hơn x). Vì vậy, bất kỳ lỗi nào trên h có thể có tác động rất lớn đến kết quả.

Vì vậy, bước đầu tiên, sửa chữa một ứng cử viên h (Tôi sẽ chỉ cho bạn trong một phút làm thế nào để chọn nó) và lấy h để tính toán số lượng h '= (x + h) - x. Nếu bạn đang sử dụng một ngôn ngữ như C, bạn phải cẩn thận để xác định h hoặc x là dễ bay hơi cho rằng tính toán không được tối ưu hóa đi.

Tiếp theo, lựa chọn h. Lỗi trong (*) có hai phần: lỗi cắt ngắn và lỗi vòng. Lỗi cắt ngắn là do công thức không chính xác:

(f(x + h) - f(x))/h = f'(x) + e1(h) 

nơi e1(h) = h/2 * sup_{x in [0,h]} |f''(x)|.

Lỗi vòng xuất phát từ thực tế là f (x + h) và f (x) gần nhau. Nó có thể được ước tính xấp xỉ như

e2(h) ~ epsilon_f |f(x)/h| 

nơi epsilon_f là độ chính xác tương đối trong tính toán của f (x) (hoặc f (x + h), đó là đóng). Điều này phải được đánh giá từ vấn đề của bạn. Đối với các chức năng đơn giản, epsilon_f có thể được lấy làm máy epsilon. Đối với những người phức tạp hơn, nó có thể tồi tệ hơn so với các đơn đặt hàng của cường độ.

Vì vậy, bạn muốn h giúp giảm thiểu e1(h) + e2(h). Cắm tất cả mọi thứ lại với nhau và tối ưu hóa trong h mang

h ~ sqrt(2 * epsilon_f * f/f'') 

mà phải được ước tính từ hàm của bạn. Bạn có thể ước tính sơ bộ. Khi nghi ngờ, hãy lấy h ~ sqrt (epsilon), nơi epsilon = độ chính xác của máy. Đối với sự lựa chọn tối ưu của h, độ chính xác tương đối mà đạo hàm được biết là sqrt (epsilon_f), tức là. một nửa số liệu quan trọng là chính xác.

Tóm lại: quá nhỏ a h => lỗi vòng, quá lớn h => lỗi cắt ngắn.

Đối với công thức bậc hai, sản lượng tính toán cùng

h ~ (6 * epsilon_f/f''')^(1/3) 

và độ chính xác thập phân của (epsilon_f)^(2/3) cho đạo hàm (mà thường là một hoặc hai chữ số có nghĩa tốt hơn so với người đầu tiên công thức đặt hàng, giả sử độ chính xác gấp đôi).

Nếu điều này quá không chính xác, vui lòng hỏi thêm phương pháp, có rất nhiều thủ thuật để có được độ chính xác cao hơn. Richardson ngoại suy là một khởi đầu tốt cho các chức năng trơn tru. Nhưng những phương pháp đó thường tính toán khá nhiều lần, điều này có thể hoặc không phải là những gì bạn muốn nếu hàm của bạn phức tạp.

Nếu bạn định sử dụng các dẫn xuất bằng số nhiều lần tại các điểm khác nhau, sẽ trở nên thú vị khi xây dựng một xấp xỉ Chebyshev.

+0

Câu trả lời rất thú vị và chi tiết, cảm ơn bạn! Tôi có thể hỏi nền của bạn là gì không? – lms

+0

@codenoob: toán học, chủ yếu. –

Các vấn đề liên quan