2

我正在尝试使用solve_ivp解决ODE系统，并且我想在每次被求解器调用时更改函数的局部变量。特别是我想更新拉格朗日乘数（lambdas_i），以便solve_ivp的下一个时间步使用前一个的值。（''reconstruct'' 函数来自一个 python 模块，该模块使用一种方法从给定时刻重建大小分布）有没有办法做到这一点？我将在下面发布代码：

import time 
import numpy as np
import scipy.integrate as integrate
from numpy import sqrt, sin, cos, pi
import math 
import pylab as pp
from pymaxent import reconstruct


start = time.time()


'''Initialize variables'''

t=np.linspace(0, 60,31)
tspan=(0,60)

initial_m=[]
for i in range(4):
    def distr(L,i=i):
        return (L**i)*0.0399*np.exp(-((L-50)**2)/200)
   
    m, err=integrate.quad(distr, 0, np.inf)
    print('m(',i,')=',m)
    initial_m.append(m)
    

''' Solving ode system using Maximum Entropy, G(L)=1+0.002*L'''

def moments(t,y):
    m0 = y[0]
    m1 = y[1]
    m2 = y[2]
    m3 = y[3]
    Lmean=m1
    σ=(m2-Lmean**2)**(1/2)
    Lmin=Lmean-3*σ
    Lmax=Lmean+4*σ
    bnds=[Lmin,Lmax]
    L=np.linspace(Lmin,Lmax)
    
    sol, lambdas_i= reconstruct(mu=y ,bnds=bnds)
    print(lambdas_i)
    
    dm0dt=0
        
    def moment1(L):
        return(sol(L)+0.002*sol(L)*L)                
    dm1dt, err1=integrate.quad(moment1,Lmin,Lmax)
    
    
    def moment2(L):
        return(2*L*(sol(L)+0.002*sol(L)*L))
    
    dm2dt, err2=integrate.quad(moment2,Lmin,Lmax)
    
    def moment3(L):
        return(3*L**2*(sol(L)+0.002*sol(L)*L))
    
    dm3dt, err3=integrate.quad(moment3,Lmin,Lmax)
    
    return(dm0dt,dm1dt,dm2dt,dm3dt)

'''Χρήση της BDF, step by step'''

r=integrate.solve_ivp(moments,tspan,initial_m,method='BDF',jac=None,t_eval=t,rtol=10**(-3))

end = time.time()

print('Total time =',{end-start})
4

1 回答 1

0

这是实现您想要的一种方法。我不会使用您的实际代码，而是使用更简单的示例问题，您可以使用相同的策略来解决您的问题。

def seq():
  l = [1, 0]
  while True:
    p = l[1]
    n = l[0] + p
    l = [p, n]
    print(n)
    input()

上面是一个示例函数，它将打印斐波那契数列的下一个（或第一个）项。（其中每一项是前两项之和。）在这种情况下，input仅用于在每次迭代之间暂停（因为它是无限的）。现在要将其转换为生成器，以提供更大的灵活性，您可以将函数重写为：

def seq():
  l = [1, 0]
  while True:
    p = l[1]
    n = l[0] + p
    l = [p, n]
    yield n

现在，如果您想获得相同的结果，您可以：

for item in seq():
  print(item)
  input()

然而，这大多是无用的。当您想要收集序列的下一个数字时，生成器的点就出现了，但在代码的任何点。不一定在循环中重复，直到你完成。你可以做到这一点：

gen = seq()
next(gen) # returns 1
next(gen) # returns 1
next(gen) # returns 2
next(gen) # returns 3
next(gen) # returns 5

等等...

解决此问题的另一种方法是使用全局变量而不是局部变量。

l = [1, 0]

def seq():
  p = l[1]
  n = l[0] + p
  l[0] = p
  l[1] = n
  return n

变量l定义在函数外部，而不是内部，所以函数退出时不会被丢弃。（以与以前相同的方式运行代码：）

while True:
  print(seq())
  input()

这两个都应该可以在您的代码中实现。

于 2021-05-17T16:16:45.213 回答