Giải phương trình vi phân phi tuyến trong python

Để tìm nghiệm của một phương trình phi tuyến tính, hãy sử dụng phương pháp chia đôi được triển khai trong tối ưu hóa mô hình con scipy. chia đôi hoặc phương pháp Newton-Raphson được triển khai trong tối ưu hóa mô hình con scipy. newton

ví dụ 3. Tìm nghiệm của phương trình phi tuyến tính

• Phương pháp chia đôi bắt đầu trên khoảng [-2, 2]___ ___ 0
• Phương pháp Newton Raphson với điểm xuất phát tại x0 = 2.

  >>> import scipy
  >>> import scipy.optimize
  >>> def f[x]:
           y = x + 2*scipy.cos[x]
           return y
           
  >>> scipy.optimize.newton[f,  2]      # starting point x0 = 2, xn+1 = xn - f[xn]/f'[xn] for n>0
  -1.0298665293222757
  

Để tìm nghiệm của một tập hợp các phương trình phi tuyến tính, tối ưu hóa mô đun con scipy. fsolve là cần thiết

Ví dụ 4. Trước tiên, hãy tìm nghiệm của phương trình phi tuyến tính một biến bằng cách sử dụng fsolve tại điểm bắt đầu  x0 = 0. 3

>>> import scipy
>>> import scipy.optimize

>>> def f[x]:
         y = x + 2*cos[x]
         return y

>>> x0 = scipy.optimize.fsolve[f, 0.3]       # 0.3 is the starting point 
>>> print x0, f[x0]
-1.0298665293

Ví dụ 5. Tìm nghiệm của hệ phương trình phi tuyến tính bắt đầu tại x0 = 1 và x1 = 1 bằng cách sử dụng fsolve

>>> import scipy
>>> import scipy.optimize

>>> def func2[x]:
        y = [x[0]*cos[x[1]] - 4, x[1]*x[0] - x[1] - 5]
        return y

>>> x0 = scipy.optimize.fsolve[func2, [1, 1]]
>>> print x0
[ 6.5041,  0.9084]

Python cung cấp một cách khác để xác định hàm bằng cách sử dụng biểu mẫu lambda. Biểu mẫu lambda cho phép tạo một đối tượng chức năng

Câu chuyện này là phần tiếp theo của câu chuyện trước đây của tôi về cách giải phương trình vi phân bằng python

Ngươi mâu

Giả sử chúng ta có tập hợp các phương trình vi phân sau

Y và P có thể là hai quần thể động vật trong cùng một khu vực. Do đó, sự tăng trưởng của một quần thể phụ thuộc vào số lượng động vật ăn thịt của các loài khác. Đối với mọi tính toán, quần thể ban đầu phải được cung cấp. Trong mô hình đầu tiên, cả hai quần thể ban đầu được đặt thành 20

Mật mã

Bây giờ chúng tôi đã xác định mô hình của mình, đã đến lúc viết chương trình python để tính toán kích thước của quần thể sau một khoảng thời gian nhất định

Đầu tiên, chúng ta cần nhập thư viện Numpy và Matplotlib. Sau đó, chúng ta có thể đặt các giá trị ban đầu. Chúng ta phải xác định một bước thời gian cho các xấp xỉ, các giá trị của Y và P ở đầu. Chúng ta cũng cần thời gian bắt đầu và kết thúc để tính số bước chúng ta phải tính

import numpy as np
import matplotlib.pyplot as pltDt = 0.01               # timestep Delta t
Y_start = 20            # initial Y
P_start = 20            # initial P
t_start = 0             # starttime
t_end = 60              # endtimen_steps = int[round[[t_end-t_start]/Dt]]    # number of timesteps

Sau đó, chúng ta có thể tạo các mảng trống để lưu trữ các giá trị được tính toán bất cứ lúc nào. Chúng ta nên điền vào phần tử 0 với các giá trị ban đầu được gán

Y_arr = np.zeros[n_steps + 1]   # create an array of zeros for Y
P_arr = np.zeros[n_steps +1]    # create an array of zeros for P
t_arr = np.zeros[n_steps + 1]   # create an array of zeros for t
t_arr[0] = t_start              # add starttime to array
Y_arr[0] = Y_start              # add initial value of Y to array
P_arr[0] = P_start              # add initial value of P to array

Bây giờ là lúc tạo một vòng lặp for sẽ lặp lại mọi dấu thời gian và tính toán các giá trị tiếp theo của P, Y và t dựa trên các phương trình vi phân của chúng. Nếu bạn muốn thử nghiệm với các phương trình vi phân khác nhau, bạn có thể làm điều đó tại đây

# Euler's method
for i in range [1, n_steps + 1]:  
   Y = Y_arr[i-1]
   P = P_arr[i-1]
   t = t_arr[i-1]
   dYdt = -0.4*Y +0.02*P*Y          # calculate the derivative of Y
   dPdt = 0.8*P - 0.01*P*P-0.1*P*Y  # calculate the derivative of Y
   Y_arr[i] = Y + Dt*dYdt  # calc. Y at next timestep,add to array
   P_arr[i] = P + Dt*dPdt  # calc. P at next timestep,add to array
   t_arr[i] = t + Dt       # add new value of t to array

Cuối cùng, chúng ta có thể vẽ biểu đồ P và Y thay đổi như thế nào theo thời gian

# plotting the result
fig = plt.figure[]                                  # create figure
plt.plot[t_arr, Y_arr, linewidth = 4, label = 'Y']    # plot Y to t 
plt.plot[t_arr, P_arr, linewidth = 4, label = 'P']    # plot P to tplt.title['Title', fontsize = 12]    # add some title to your plot
plt.xlabel['t [in seconds]', fontsize = 12]
plt.ylabel['Y[t], P[t]', fontsize = 12]
plt.xticks[fontsize = 12]
plt.yticks[fontsize = 12]
plt.grid[True]                        # show grid
plt.axis[[t_start, t_end, 0, 50]]     # show axes measures
plt.legend[]
plt.show[]

Như bạn có thể thấy, một trạng thái cân bằng đã đạt được. Bây giờ chúng ta đã giải hệ phương trình vi phân này. Bạn có thể chơi xung quanh với mã này và thay đổi một số biến và xem kết quả thay đổi như thế nào hoặc nếu đạt đến điểm cân bằng hoặc thậm chí có thể là một trạng thái cân bằng khác. Bạn vẫn nên đảm bảo chọn bước thời gian dt thích hợp

giản đồ pha

Nếu chúng ta vẽ biểu đồ Y theo P thay vì theo t, chúng ta sẽ có được cái gọi là giản đồ pha, khá hữu ích cho việc nghiên cứu các trạng thái cân bằng. Trong trường hợp này, nó không thú vị bằng, nhưng có thể tạo ra một hình ảnh đẹp. Trong trường hợp này, cả ba kết hợp ban đầu hội tụ tại một điểm duy nhất. Cái này được gọi là bồn rửa xoắn ốc

Biểu đồ pha cho hệ phương trình vi phân với các giá trị ban đầu trong truyền thuyết

Nếu bạn đã hiểu mã này và các lý thuyết hỗ trợ nó, bạn sẽ có cơ sở tuyệt vời để giải bất kỳ hệ phương trình vi phân nào bằng phương pháp số

Bạn có thể giải phương trình vi phân bằng Python không?

Có thể giải phương trình vi phân phi tuyến tính không?

Phương pháp nào là tốt nhất để giải phương trình phi tuyến tính?

Sự khác biệt giữa Odeint và giải quyết_ivp là gì?