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 -2sage: 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_rk4sage: 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)
洛伦兹方程：
令 σ = 10 sigma=10 σ=10, ρ = 28 rho=28 ρ=28, β = 3 beta=3 β=3
0时刻初值给的(1,1,1)
sage: from sage.calculus.desolvers import desolve_system_rk4sage: x,y,z,t=var('x y z t')sage: p=10sage: r=28sage: b=3sage: 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中有所描述，下面的代码也是参考了文档中的交互式代码改过来的，可以直接运行出来。
令 μ = 2 mu=2 μ=2, σ = 0.5 sigma=0.5 σ=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 = 2sigma = 0.5plots_at_a_time = 10number_of_steps = 99clear_plot = Trueemplot = 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中关于随机微分方程的部分