python中的数值延续
pacop的Python项目详细描述
压缩
pacopy为python中的pde提供了各种numerical parameter continuation算法。
pacopy是后端不可知论者,所以你的问题是否由
SciPy,FEniCS,
pyfvm,或任何其他python包。唯一的事
用户必须提供一个具有一些简单方法的类,例如函数求值
f(u, lmbda)
,雅可比矩阵a解算器jacobian_solver(u, lmbda, rhs)
等。
有些pacopy文档是可用的here。
示例
布拉图经典的Bratu problem 一维dirichlet边界条件。要再现情节,你首先必须 指定问题;这是经典的有限差分近似:
classBratu1d(object):def__init__(self):self.n=51h=1.0/(self.n-1)self.H=numpy.full(self.n,h)self.H[0]=h/2self.H[-1]=h/2self.A=(scipy.sparse.diags([-1.0,2.0,-1.0],[-1,0,1],shape=(self.n,self.n))/h**2)returndefinner(self,a,b):"""The inner product of the problem. Can be numpy.dot(a, b), but factoring in the mesh width stays true to the PDE. """returnnumpy.dot(a,self.H*b)defnorm2_r(self,a):"""The norm in the range space; used to determine if a solution has been found. """returnnumpy.dot(a,a)deff(self,u,lmbda):"""The evaluation of the function to be solved """out=self.A.dot(u)-lmbda*numpy.exp(u)out[0]=u[0]out[-1]=u[-1]returnoutdefdf_dlmbda(self,u,lmbda):"""The function's derivative with respect to the parameter. Used in Euler-Newton continuation. """out=-numpy.exp(u)out[0]=0.0out[-1]=0.0returnoutdefjacobian_solver(self,u,lmbda,rhs):"""A solver for the Jacobian problem. """M=self.A.copy()d=M.diagonal().copy()d-=lmbda*numpy.exp(u)M.setdiag(d)# Dirichlet conditionsM.data[0][self.n-2]=0.0M.data[1][0]=1.0M.data[1][self.n-1]=1.0M.data[2][1]=0.0returnscipy.sparse.linalg.spsolve(M.tocsr(),rhs)
然后将对象传递到Pacopy的任何方法,例如Euler-Newton(ArcLength) 续:
problem=Bratu1d()# Initial guessu0=numpy.zeros(problem.n)# Initial parameter valuelmbda0=0.0lmbda_list=[]values_list=[]defcallback(k,lmbda,sol):# Use the callback for plotting, writing data to files etc.fig=plt.figure()ax1=fig.add_subplot(111)ax1.set_xlabel("$\\lambda$")ax1.set_ylabel("$||u||_2$")ax1.grid()lmbda_list.append(lmbda)values_list.append(numpy.sqrt(problem.inner(sol,sol)))ax1.plot(lmbda_list,values_list,"-x",color="#1f77f4")ax1.set_xlim(0.0,4.0)ax1.set_ylim(0.0,6.0)return# Natural parameter continuation# pacopy.natural(problem, u0, lmbda0, callback, max_steps=100)pacopy.euler_newton(problem,u0,lmbda0,callback,max_steps=500,newton_tol=1.0e-10)
金茨堡-兰道
Ginzburg-Landau equations模型 极端ii型超导体在磁场下的行为。上面的例子 (详细内容见 here中) 显示磁场强度的参数连续性。上的情节 右侧显示了使用 cplot。
安装
pacopy是available from the Python Package Index,所以只需键入
pip install -U pacopy
安装或升级。
测试
要运行pacopy单元测试,请签出此存储库并键入
pytest
分布
创建新版本
把
__version__
号撞一下,发布到pypi和github:
make publish
许可证
pacopy在MIT license下发布。