欢迎您访问365答案网,请分享给你的朋友!
生活常识 学习资料

SageMath的使用

时间:2023-05-24
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 -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中关于随机微分方程的部分

Copyright © 2016-2020 www.365daan.com All Rights Reserved. 365答案网 版权所有 备案号:

部分内容来自互联网,版权归原作者所有,如有冒犯请联系我们,我们将在三个工作时内妥善处理。