Tôi phải giải các phương trình siêu việt sau cos(x)/x=c cho không đổi c. Ví dụ: tôi đã thực hiện một mã ngắn trong Mathicala, trong đó tôi đã tạo một danh sách các giá trị ngẫu nhiên cho hằng số C I did a short code in Mathematica, where I generated a list of random values for constant c const = Table[RandomReal[{0, 5}], {i, 1, 10}]
(*{1.67826, 0.616656, 0.290878, 1.10592, 0.0645222, 0.333932, 3.59584, \
2.70337, 3.91535, 2.78268}*)
Hơn tôi đã xác định chức năng f[x_, i_] := Cos[x]/x - const[[i]]
và bắt đầu tìm kiếm rễ: Table[FindRoot[f[x, i] == 0, {x, 0.1}][[1, 2]], {i, 1, Length[const]}]
(*{0.517757, 0.947103, 1.21086, 0.694679, 1.47545, 1.16956, 0.26816, \
0.347764, 0.247615, 0.338922}*)
Bây giờ tôi rất thích lập trình một cái gì đó tương tự trong Python (có thể sử dụng Numpy?) Nhưng tôi thực sự không thể tìm thấy bất kỳ câu trả lời nào tốt cho một vấn đề như thế. Ai đó có thể giúp đỡ? Phương trình siêu việt mà chúng ta sẽ giải quyết trong ví dụ này là: $ \ ln \ left (\ dfrac {\ dot {q} rt} {\ dot {v} p_0l_ \ mathrm {m}} \ right) = \ dfrac {l_ \ mathrm {m}} dfrac {1} {t_0}-\ dfrac {1} {t} \ right) $ Phương trình này xác định nhiệt độ cân bằng của chất lỏng sôi đang được làm mát bằng cách làm mát bay hơi. Áp suất hơi của nó giảm thông qua một máy bơm với tốc độ bơm $ \ dot {v} $ và chất lỏng đang hấp thụ tải nhiệt $ \ dot {q} $ từ môi trường xung quanh. Sự phụ thuộc của nhiệt độ sôi vào áp suất hơi được xác định từ mối quan hệ calusius-Clapeyron. Chất lỏng có nhiệt độ cân bằng $ T_0 $ khi áp suất hơi là $ P_0 $. $ R $ là hằng số khí phổ quát $ l_ \ mathrm {m} $ là nhiệt tiềm ẩn mol
$L_\mathrm{m}$ is the molar latent heat
- Nhập một số giá trị số của các đại lượng khác nhau. Mục tiêu là tìm nhiệt độ cân bằng không thể giải quyết được để phân tích.
Nhiệt độ cân bằng cho tải nhiệt $ \ dot {q} = 2.5 ~ \ mathrm {w} $ và tốc độ bơm là $ \ dot {v} = 1.0 \ lần 10^{-3} ~ \ mathrm {m} ^2/\ mathrm {s} $ là $ t = 1.291 ~ \ mathrm {k} $. Chúng ta có thể xem xét nếu giải pháp này có ý nghĩa. Chúng tôi thực hiện việc này âm mưu của tôi ở bên trái và bên phải của phương trình siêu việt như là một hàm của $ t $ trên cùng một biểu đồ và tìm kiếm giao điểm. - Đầu tiên xác định các hàm cho LHS & RHS của phương trình siêu việt.
Vì đây là một phương trình siêu việt, việc tìm kiếm tất cả rễ không phải là một lựa chọn. Bạn không thể (nói chung) tìm một biểu thức đưa ra dạng đóng cho rễ của phương trình. Vì vậy, bạn cần phải làm một số phân tích trước. Vẽ vẽ hàm, hoặc vẽ sơ đồ phía bên trái và bên phải của phương trình (và do đó, các giao điểm của cả hai đường cong là rễ) sẽ cho bạn một số gợi ý về cách giải quyết vấn đề. Tôi định nghĩa các LHS là $ \ tan (x) $ và rhs là $ \ frac {2x} {x^{2} -1} $. Những gì bạn có thể quan sát là như sau: - Vì cả LHS và RHS đều đối xứng (chống) xung quanh $ x = 0 $, bạn chỉ cần xem xét miền $ x $ dương.
- Không có gốc trong khoảng $ [0,1 [$ ($ x = 1 $ là vị trí của cực của hàm RHS).
- Có một gốc trong khoảng $ [1, \ frac {\ pi} {2}] $.
- Vì chức năng của RHS giảm một cách đơn điệu với $ x> 1 $, sẽ có một gốc duy nhất trong mỗi khoảng $ [(2K-1) \ frac {\ pi} {2}, (2K+1) \ frac { \ pi} {2}] $ cho $ k = 1,2,3, \ ldots $.
- Bạn có thể thấy rằng hàm RHS chuyển đến 0 với giá $ x \ đến \ infy $. Điều này ngụ ý rằng với $ x $ lớn hơn, các giải pháp sẽ hội tụ với các giao cắt bằng 0 của hàm tiếp tuyến (tức là $ r_ {l} \ xấp xỉ (l+1) \ pi $, trong đó $ r_ {l} $ là $ l $ '-r gốc. Bạn có thể sử dụng các ước tính này như dự đoán ban đầu cho một phương thức lặp.
Vì vậy, những gì sẽ là chiến lược để tìm thấy (không phải tất cả, mà là một số lượng tùy ý) của rễ? Bạn sử dụng một phương pháp giải quyết lặp trong mỗi khoảng. Bây giờ bạn có sự lựa chọn giữa rất nhiều phương pháp nhưng hai phương pháp có vẻ rất phù hợp: - Vì bạn có một ước tính rất tốt cho gốc (một phỏng đoán ban đầu tốt), bạn có thể sử dụng Newton-Raphson. Nhưng sau đó bạn cần tính toán đạo hàm (trong trường hợp này, bạn có thể làm điều đó một cách phân tích).
- Vì mỗi gốc được biết là trong một khoảng thời gian, bạn có thể sử dụng phương pháp không có đạo hàm của Dekker-Brent. Chỉ cần đảm bảo rằng bạn không sử dụng các ranh giới nghiêm ngặt cho mỗi khoảng thời gian, bởi vì chúng là cực của hàm tiếp tuyến và Dekker-Brent sẽ không biết cách xử lý điều đó. Thêm một số lề vào khoảng thời gian ban đầu.
Sử dụng mã Python bên dưới, tôi có thể tìm thấy rễ $ n $ đầu tiên bằng cách sử dụng phương pháp này. Một âm mưu cho khoảng [200.300] được đưa ra làm ví dụ. Vì tôi lười biếng và không muốn tính toán đạo hàm, tôi đã chọn tùy chọn Dekker-Brent (nơi tôi yêu cầu dung sai tương đối là $ 1E-14 $. Trên máy tính xách tay Intel i5. Tất nhiên, bạn có thể xác nhận kết quả bằng cách tính toán dư lượng (giá trị hàm trong mỗi gốc được tính toán). Mã Python: import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import brentq
plt.figure(figsize=(12,8))
plt.style.use('bmh')
x = np.linspace(0,20,1002)
plt.ylim(-5,5)
plt.plot(x,np.tan(x),label='$tan(x)$')
plt.plot(x,2*x/(x*x-1),label='$2x/(x^2-1)$')
plt.legend(loc=1)
plt.savefig('fun.png',dpi=300,bbox_inches = 'tight')
def f(x):
# The function
return np.tan(x)-2*x/(x*x-1)
def roots(N):
# Find N roots of the equation
#Allocate space
roots = np.zeros(N)
#Margin to stay away from poles
margin = 1e-8
#First root
roots[0] = brentq(f, 1.0 + margin, np.pi/2 - margin, rtol=1e-14)
#Subsequent N-1 roots
for i in range(1,N):
left = (2*i - 1)*np.pi/2
right = (2*i + 1)*np.pi/2
roots[i] = brentq(f, left + margin, right - margin, rtol = 1e-14)
return roots
result = roots(1000)
left = 200
right = 300
plt.figure(figsize=(24,8))
plt.style.use('bmh')
x = np.linspace(left,right,(right-left)*101)
plt.ylim(0,2e-2)
plt.plot(x,np.tan(x),label='$tan(x)$')
plt.plot(x,2*x/(x*x-1),label='$2x/(x^2-1)$')
rts = result[result>=left]
rts = rts[rts<=right]
plt.plot(rts,np.tan(rts),'o', markersize=14, fillstyle='none', label='roots')
plt.legend(loc=1)
plt.savefig('roots1.png',dpi=300,bbox_inches = 'tight')
| 
