分數階累加的Python實現

分數階累加是分數階差分的逆運算，它不僅可用于分數階差分方程的分析 ，也可以用于建立分數階灰色模型。然而許多初學者在動手實現分數階灰色模型時經常發現非常困難，究其原因其實是對定義公式的分析不夠，對相應程序語言的特性不熟悉。

本文將從分數階累加的定義出發，深入分析其計算過程，結合Python語言的特性，詳細講解其實現過程。

1、 分數階累加的定義

對任意原始序列 $x^{(0)}(1),x^{(0)}(2),...,x^{(0)}(n)$，其分數階累加定義為：
$$x^{(r)}(k)=\sum _{i=1}^k\left(\begin{array}{c}k?i+r?1\\ k?i\end{array}\right)x^{(0)}(i)\text{\hspace{1em}(1)}$$

其中：
$$\left(\begin{array}{c}k?i+r?1\\ k?i\end{array}\right)=\frac{(r+k?i?1)(r+k?i?2)?(r+1)r}{(k?i)!}\text{\hspace{1em}(2)}$$

通俗地講，分數階累加實際上就是一種加權求和，而在具體實現時關鍵是求出其加權系數 (2)，該系數稱為廣義二項式系數或廣義牛頓二項式系數。通常初學者在看到該定義的第一反應即是用循環來實現，這種方式直接簡單，但效率較低。尤其在使用類似Python這類腳本語言時其速度非常之慢。

2、廣義牛頓二項式系數的計算

要實現分數階累加的快速算法，關鍵的問題就在于如何快速算出其加權系數。對于一些不常用科學計算類程序的人而言，通常不太容易想到簡單的解決方法。事實上，Python中提供了一個函數binom正好可以解決該問題。

binom函數是Scipy庫中的一個函數。在Scipy的官方文檔中對該函數并沒有過多的介紹，那么我們先來簡單寫幾個測試看看它的用法：

from scipy.special import binom# integers print(binom(10,3))# decimals print(binom(10,0.3))# array_like data print(binom([5,8,10],[0.1,2,3.1]))

運行結果：

120.0 2.246507496036618 [ 1.24554362 28. 129.20102313]

從上述代碼中不難看出，binom函數支持的輸入是十分簡單的，第一個變量是二項式系數中的n，第二個變量是k，即是在n中選出k個的組合數。同時它也支持小數和數組，這為我們進一步簡化計算過程提供了極大的便利。

3、分數階累加的矩陣形式及向量計算

本文既然是要討論快速算法，那么自然就會盡量避免循環。雖然在目前看來，完全可以直接使用binom函數對其進行計算，但這樣仍然要使用2層循環，并且計算$\frac{n(n+1)}{2}$次。試想假設我們有10個點要計算，那么就得調用5050次，這個量級對于Python的循環是很低效的。

實際上Python中的numpy庫是具有十分完善的矩陣計算功能，其計算效率非常高。我們要想辦法利用這一特性，那么很自然地就需要引入它的矩陣定義：

$$\left[\begin{array}{c}{x}^{(r)}(1),x^{(r)}(2),?,x^{(r)}(n)\end{array}\right]=\left[\begin{array}{c}{x}^{(0)}(1),x^{(0)}(2),?,x^{(0)}(n)\end{array}\right]?\begin{array}{ccccc}1& \left(\begin{array}{c}r\\ 1\end{array}\right)& ?& \left(\begin{array}{c}r+n?3\\ n?2\end{array}\right)& \left(\begin{array}{c}r+n?2\\ n?1\end{array}\right)\\ 0& 1& ?& \left(\begin{array}{c}r+n?4\\ n?3\end{array}\right)& \left(\begin{array}{c}r+n?3\\ n?2\end{array}\right)\\ ?& ?& ?& ?& ?\\ 0& 0& ?& 1& \left(\begin{array}{c}r\\ 1\end{array}\right)\\ 0& 0& ?& 0& 1\end{array}?\text{\hspace{1em}(3)}$$

這里有一個特別要注意的點：處在同一右對角線上的所有元素都相等。進一步我們可以發現一個更直接有利于我們編程的事實：從右向左的每一列都相當于其右列去掉最后一個元素，再向左移動一位。那么我們容易得出一個非常簡單且有用的結論：我們在編程實現時只需要計算最后一列。 此時的時間復雜度已經從$O(n^2)$降到了線性時間復雜度$O(n)$，同時，其內存消耗也由原來的$n^2$單位變成了$n$個單位。

要計算最后一列系數，實際上只需要考慮定義式 (3) 中$k=n$的情形，利用Python的向量計算格式，可以很容易實現這一目的：

import numpy as np from scipy.special import binomr = 0.1 k = 10 i = np.linspace(1,k,k)fnum = k - i + r -1 ki = k - icoefs = binom(fnum,ki)print(coefs)

運行結果：

[0.01447564 0.01608405 0.01812287 0.02079674 0.02446675 0.02983750.0385 0.055 0.1 1. ]

4、分數階累加的實現

結合上述分析，這里我們得到的系數coefs實際上是矩陣式 (3) 中的第一行。再次觀察矩陣式注意到在計算累加時是采用原序列的行左乘矩陣的列的方式進行。從編程的角度考慮則不難發現要計算第$k$個分數階累加值，即是將前$k$個原序列點中的值組成的向量對應乘以系數的倒數前$k$個值再求和。這一步實現的時候就需要用到一層循環了。

x = np.random.rand(10) xr = np.zeros(10) for k in range(10):xr[k] = sum(x[:k+1]*coefs[-(k+1):])print(x) print(xr)

運行結果：

[0.42343781 0.53237537 0.75292602 0.32316064 0.6931979 0.496077060.74748506 0.3630722 0.29325613 0.95218414] [0.42343781 0.57471915 0.82945264 0.44403624 0.80005567 0.638403230.89195739 0.5386026 0.45048115 1.09707706]

5、分數階累加的完整代碼

最后為方便使用，我們將該方法封裝為完整的函數

from scipy.special import binom import numpy as np# compute the r-order accumulation def fago(x,r):n = len(x)i = np.linspace(1,n,n)coefs = binom(n - i + r -1,n - i)xr = np.zeros(10)for k in range(10):xr[k] = sum(x[:k+1]*coefs[-(k+1):])return xr

簡單測試一下這個函數：

import time t1 = time.time() print(fago(x,0.1)) t2 = time.time()print('Runtime:',t2-t1,'s')

運行結果：

array([0.42343781, 0.57471915, 0.82945264, 0.44403624, 0.80005567,0.63840323, 0.89195739, 0.5386026 , 0.45048115, 1.09707706]) Runtime: 0.0010225772857666016 s

可以看到和剛才的結果一模一樣，運行時間只有0.001秒。

更多關于分數階累加的應用可查看相應論文：

A_novel_fractional_grey_system_model_and_its_application

A novel fractional time delayed grey model with Grey Wolf Optimizer and its applications in forecasting the natural gas and coal consumption in Chongqing China

The novel fractional discrete multivariate grey system model and its applications

總結

以上是生活随笔為你收集整理的分数阶累加的Python实现的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。