我如何简化代码以获取Jacobian Pass

发布于 2025-02-03 03:51:24 字数 2673 浏览 2 评论 0原文

我正在尝试计算一个5x12矩阵,该矩阵包含相对于12个参数的5个函数的衍生物(S,D,Beta,...)。我正在使用Jacobian来计算事物,但是我在此处遇到此错误:()缺少1所需的位置参数:'x'。非常喜欢一些指针。这是代码

import autograd.numpy as np
import matplotlib.pyplot as plt
 
def HCVModel(x,params):
    s=params["s"]
    d=params["d"]
    beta=params["beta"]
    a=params["a"]
    p=params["p"]
    k=params["k"]
    mu=params["mu"]
    q=params["q"]
    g=params["g"]
    c=params["c"]
    h=params["h"]
    b=params["b"]
 
    xdot=np.array([s-d*x[0]-beta*x[2]*x[0],beta*x[2]*x[0]-a*x[1]-p*x[1]*x[3],k*x[1]-mu*x[2]-q*x[2]*x[4],c*x[1]*x[3]-b*x[3],g*x[2]*x[4]-h*x[4]])
 
    return xdot
 
def RK4(f,x0,t0,tf,dt):
    t=np.arange(t0,tf,dt)
    nt=t.size
    nx=x0.size
    x=np.zeros((nx,nt))
    x[:,0]=x0
    for k in range(nt-1):
        k1=dt*f(t[k],x[:,k])
        k2=dt*f(t[k]+dt/2,x[:,k]+k1/2)
        k3=dt*f(t[k]+dt/2,x[:,k]+k2/2)
        k4=dt*f(t[k]+dt,x[:,k]+k3)
        dx=(k1+2*k2+2*k3+k4)/6
        x[:,k+1]=x[:,k]+dx
    return x,t
 
f = lambda t, x : HCVModel(x,params)
params={"s":10*10**4,"d":10*10**-1,"beta":0,"a":10*10**-3,"p":10*10**-3,"k":10*10**-3,"mu":10*10**-3,"q":10*10**-3,"g":0,"c":0,"h":10*10**-3,"b":10*10**-3}
x0=np.array([10000,2,3,1,1.5])
t0=0
tf=100
dt=10**-2
 
x, t = RK4(f,x0,t0,tf,dt)
 
 
 
from autograd.scipy.integrate import odeint
from autograd import jacobian
 
s=10*10**4
d=10*10**-1
beta=0
a=10*10**-3
p=10*10**-3
k=10*10**-3
mu=10*10**-3
q=10*10**-3
g=0
c=0
h=10*10**-3
b=10*10**-3
 
dCdteta=jacobian(f,0)
teta_sensitivity=dCdteta(np.array([s,d,beta,a,p,k,mu,q,g,c,h,b]))
s_sensitivity=teta_sensitivity[:0,0]
d_sensitivity=teta_sensitivity[:0,1]
beta_sensitivity=teta_sensitivity[:0,2]
a_sensitivity=teta_sensitivity[:0,3]
p_sensitivity=teta_sensitivity[:0,4]
k_sensitivity=teta_sensitivity[:0,5]
mu_sensitivity=teta_sensitivity[:0,6]
q_sensitivity=teta_sensitivity[:0,7]
g_sensitivity=teta_sensitivity[:0,8]
c_sensitivity=teta_sensitivity[:0,9]
h_sensitivity=teta_sensitivity[:0,10]
b_sensitivity=teta_sensitivity[:0,11]
plt.plot(tspan, np.abs(fds), label='fd dC/ds')
plt.plot(tspan, np.abs(fdd), label='fd dC/dd')
plt.plot(tspan, np.abs(fdd), label='fd dC/dd')
plt.plot(tspan, np.abs(fdbeta), label='fd dC/dbeta')
plt.plot(tspan, np.abs(fda), label='fd dC/da')
plt.plot(tspan, np.abs(fdp), label='fd dC/dp')
plt.plot(tspan, np.abs(fdk), label='fd dC/dk')
plt.plot(tspan, np.abs(fdmu), label='fd dC/dmu')
plt.plot(tspan, np.abs(fdq), label='fd dC/dq')
plt.plot(tspan, np.abs(fdg), label='fd dC/dg')
plt.plot(tspan, np.abs(fdc), label='fd dC/dc')
plt.plot(tspan, np.abs(fdh), label='fd dC/dh')
plt.plot(tspan, np.abs(fdb), label='fd dC/db')
plt.xlabel('t')
plt.ylabel('sensitivity')
plt.show()

I'm trying to compute a 5x12 matrix containing the derivatives of 5 functions with respect to 12 parameters (s,d,beta,...). I'm using the jacobian to compute the thing but I get this error here TypeError: () missing 1 required positional argument: 'x'. Would very much like some pointers. Here is the code

import autograd.numpy as np
import matplotlib.pyplot as plt
 
def HCVModel(x,params):
    s=params["s"]
    d=params["d"]
    beta=params["beta"]
    a=params["a"]
    p=params["p"]
    k=params["k"]
    mu=params["mu"]
    q=params["q"]
    g=params["g"]
    c=params["c"]
    h=params["h"]
    b=params["b"]
 
    xdot=np.array([s-d*x[0]-beta*x[2]*x[0],beta*x[2]*x[0]-a*x[1]-p*x[1]*x[3],k*x[1]-mu*x[2]-q*x[2]*x[4],c*x[1]*x[3]-b*x[3],g*x[2]*x[4]-h*x[4]])
 
    return xdot
 
def RK4(f,x0,t0,tf,dt):
    t=np.arange(t0,tf,dt)
    nt=t.size
    nx=x0.size
    x=np.zeros((nx,nt))
    x[:,0]=x0
    for k in range(nt-1):
        k1=dt*f(t[k],x[:,k])
        k2=dt*f(t[k]+dt/2,x[:,k]+k1/2)
        k3=dt*f(t[k]+dt/2,x[:,k]+k2/2)
        k4=dt*f(t[k]+dt,x[:,k]+k3)
        dx=(k1+2*k2+2*k3+k4)/6
        x[:,k+1]=x[:,k]+dx
    return x,t
 
f = lambda t, x : HCVModel(x,params)
params={"s":10*10**4,"d":10*10**-1,"beta":0,"a":10*10**-3,"p":10*10**-3,"k":10*10**-3,"mu":10*10**-3,"q":10*10**-3,"g":0,"c":0,"h":10*10**-3,"b":10*10**-3}
x0=np.array([10000,2,3,1,1.5])
t0=0
tf=100
dt=10**-2
 
x, t = RK4(f,x0,t0,tf,dt)
 
 
 
from autograd.scipy.integrate import odeint
from autograd import jacobian
 
s=10*10**4
d=10*10**-1
beta=0
a=10*10**-3
p=10*10**-3
k=10*10**-3
mu=10*10**-3
q=10*10**-3
g=0
c=0
h=10*10**-3
b=10*10**-3
 
dCdteta=jacobian(f,0)
teta_sensitivity=dCdteta(np.array([s,d,beta,a,p,k,mu,q,g,c,h,b]))
s_sensitivity=teta_sensitivity[:0,0]
d_sensitivity=teta_sensitivity[:0,1]
beta_sensitivity=teta_sensitivity[:0,2]
a_sensitivity=teta_sensitivity[:0,3]
p_sensitivity=teta_sensitivity[:0,4]
k_sensitivity=teta_sensitivity[:0,5]
mu_sensitivity=teta_sensitivity[:0,6]
q_sensitivity=teta_sensitivity[:0,7]
g_sensitivity=teta_sensitivity[:0,8]
c_sensitivity=teta_sensitivity[:0,9]
h_sensitivity=teta_sensitivity[:0,10]
b_sensitivity=teta_sensitivity[:0,11]
plt.plot(tspan, np.abs(fds), label='fd dC/ds')
plt.plot(tspan, np.abs(fdd), label='fd dC/dd')
plt.plot(tspan, np.abs(fdd), label='fd dC/dd')
plt.plot(tspan, np.abs(fdbeta), label='fd dC/dbeta')
plt.plot(tspan, np.abs(fda), label='fd dC/da')
plt.plot(tspan, np.abs(fdp), label='fd dC/dp')
plt.plot(tspan, np.abs(fdk), label='fd dC/dk')
plt.plot(tspan, np.abs(fdmu), label='fd dC/dmu')
plt.plot(tspan, np.abs(fdq), label='fd dC/dq')
plt.plot(tspan, np.abs(fdg), label='fd dC/dg')
plt.plot(tspan, np.abs(fdc), label='fd dC/dc')
plt.plot(tspan, np.abs(fdh), label='fd dC/dh')
plt.plot(tspan, np.abs(fdb), label='fd dC/db')
plt.xlabel('t')
plt.ylabel('sensitivity')
plt.show()

如果你对这篇内容有疑问,欢迎到本站社区发帖提问 参与讨论,获取更多帮助,或者扫码二维码加入 Web 技术交流群。

扫码二维码加入Web技术交流群

发布评论

需要 登录 才能够评论, 你可以免费 注册 一个本站的账号。
列表为空,暂无数据
我们使用 Cookies 和其他技术来定制您的体验包括您的登录状态等。通过阅读我们的 隐私政策 了解更多相关信息。 单击 接受 或继续使用网站,即表示您同意使用 Cookies 和您的相关数据。
原文