# -*- coding: utf-8 -*- ''' Time: XXXX author: XX Project: XXXX File: XXX '''import numpy as np from matplotlib import pyplot as plt import seaborn as sns#由題目定義矩陣 A = np.matrix([[0, 1], [0, 0]]) AT = np.matrix([[0, 0], [1, 0]]) B = np.matrix([[0], [1]]) BT = np.matrix([[0, 1]]) F = np.matrix([[1, 0], [0, 2]]) Q = np.matrix([[2, 1], [1, 4]]) R = np.matrix([[1/2]]) RN = 2#根據(jù)邊界條件給定的值 step_num = 100 t = 3 step = -t / step_num P = F# 定義黎卡提方程 def Ricatti_P(t, P):f = -(P * A + A.T * P - P * B * R.I * B.T * P + Q)return fys_0, ys_1, ys_2, ys_3 = [], [], [], [] ts = [] while t > 0:t += stepk1 = step * Ricatti_P(t, P)k2 = step * Ricatti_P(t + step * 0.5, P + k1 * step * 0.5)k3 = step * Ricatti_P(t + step * 0.5, P + k2 * step * 0.5)k4 = step * Ricatti_P(t + step, P + k3 * step)P = P + (k1 + k2 * 2 + k3 * 2 + k4)/6P = np.array(P)ts.append(t)ys_0.append(P[0][0])ys_1.append(P[0][1])ys_2.append(P[1][0])ys_3.append(P[1][1])print(ys_0)print(P)sns.set() fig, axes = plt.subplots(2, 2) axes[0][0].plot(ts, ys_0, label='p1') axes[0][1].plot(ts, ys_1, label='p2') axes[1][0].plot(ts, ys_2, label='p3') axes[1][1].plot(ts, ys_3, label='p4') for i in range(2):for j in range(2):axes[i][j].set_xlabel('t')axes[i][j].legend(loc=3) plt.suptitle('Riccati equation\'s solve ', size=18) plt.show()

總結(jié)

以上是生活随笔為你收集整理的Python实现四阶龙格库塔法求解Ricatti方程的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。