SageMath的使用，包括矩阵、微分方程和随机过程

sage的官方网站是https://www.sagemath.org/
Sage和Python紧密相连，正好一并学习。

矩阵特征值求解

利用sage可以生成随机矩阵：

sage: m = random_matrix(RDF,10)
print(m)

生成出随机矩阵，N=10
进而可以利用其eigenvalues来求本征值：

sage: e = m.eigenvalues()
sage: w = [(i, abs(e[i])) for i in range(len(e))]
sage: show(points(w))

sage对于大矩阵的求解比较方便
大矩阵：

sage: m = random_matrix(RDF,500)
sage: e = m.eigenvalues()  
sage: w = [(i, abs(e[i])) for i in range(len(e))]
sage: show(points(w))

根据上面的代码，可以求出特征值的分布。这里N=500。

---

用sage求常微分方程和方程组

一个普通的常微分

sage: y = function('y')(x)
sage: de = diff(y,x) + y -2
sage: h = desolve_rk4(de, y, step=.05, ics=[0,3])

sage: h1 = desolve(de, y, ics=[0,3])
sage: plot(h1,(x,0,5),color='red')+points(h)

下面试图常微分方程组：

from sage.calculus.desolvers import desolve_system_rk4
sage: x,y,t=var('x y t')
sage: P=desolve_system_rk4([x*(1-y),-y*(1-x)],[x,y],ics=[0,0.5,2],ivar=t,end_points=20)
sage: Q=[ [i,j] for i,j,k in P]
sage: LP=list_plot(Q)

sage: Q=[ [j,k] for i,j,k in P]
sage: LP=list_plot(Q)

points(P)

洛伦兹方程：

令$\sigma =10$,$\rho =28$,$\beta =3$
0时刻初值给的(1,1,1)

sage: from sage.calculus.desolvers import desolve_system_rk4
sage: x,y,z,t=var('x y z t')
sage: p=10
sage: r=28
sage: b=3
sage: P=desolve_system_rk4([p*(y-x), x*(r-z)-y, x*y-b*z],[x,y,z],ics=[0,1,1,1],ivar=t,end_points=20,step=0.01)
sage: print(P) 

points(P)

这里用的是4阶龙格库塔法进行的数值求解

随机过程（布朗运动）

描述布朗运动的微分方程在https://wiki.sagemath.org/interact/diffeq中有所描述，下面的代码也是参考了文档中的交互式代码改过来的，可以直接运行出来。

令$\mu =2$,$\sigma =0.5$，
可以在sage中用欧拉-丸山法求随机微分方程：

def EulerMaruyama(xstart, ystart, xfinish, nsteps, f1, f2):
    sol = [ystart]
    xvals = [xstart]
    h = N((xfinish-xstart)/nsteps)
    for step in range(nsteps):
        sol.append(sol[-1] + h*f1(sol[-1]) + h^(.5)*f2(sol[-1])*normalvariate(0,1)) 
        xvals.append(xvals[-1] + h)
    return list(zip(xvals,sol))

out = Graphics()
save(out,'temp')

mu = 2
sigma = 0.5
plots_at_a_time = 10
number_of_steps = 99
clear_plot = True
emplot = list_plot(EulerMaruyama(0,1,1,number_of_steps,lambda x: mu*x,lambda x:sigma*x),plotjoined=True)
for i in range(1,plots_at_a_time):
    emplot = emplot + list_plot(EulerMaruyama(0,1,1,100,lambda x: mu*x,lambda x:sigma*x),plotjoined=True)
if clear_plot:
    out = emplot
    save(out,'temp')
else:
    out = load('temp')
    out = out + emplot
    save(out,'temp')
show(out, figsize = [8,5])

参考的文档是https://wiki.sagemath.org/interact/diffeq中关于随机微分方程的部分