在python中实现Backward Euler方法来解决摆锤

Question

我正在尝试为摆F和dF设置一个隐式求解器，其定义为：

def func(y,t):
    ### Simplified the Function to remove friction since it canceled out

    x,v = y[:3],y[3:6]
    grav = np.array([0.,  0., -9.8 ])
    lambd = (grav.dot(x)+v.dot(v))/x.dot(x)

    return np.concatenate([v, grav - lambd*x] )
def dF_matrix(y):

    n=y.size
    dF=np.zeros((6,6))

    xp=np.array([y[3],y[4],y[5]])[np.newaxis]

    mass=1.
    F1=0.
    F2=0.
    F3=-mass*9.8
    F=np.array([F1,F2,F3])[np.newaxis]

    phix=2.*y[0]
    phiy=2.*y[1]
    phiz=2.*y[2]
    G=np.array([phix,phiy,phiz])[np.newaxis]

    H=2.*np.eye(3)
    lambd=(mass*np.dot(xp,np.dot(H,xp.T))+np.dot(F,G.T))/np.dot(G,G.T)

    dF[0,3]=1
    dF[1,4]=1
    dF[2,5]=1
    dF[3,0]=(y[0]*F1+2*lambd)/mass
    dF[3,1]=(y[0]*F2)/mass
    dF[3,2]=(y[0]*F3)/mass
    dF[3,3]=phix*y[3]
    dF[3,4]=phix*y[4]
    dF[3,5]=phix*y[5]
    dF[4,0]=(y[1]*F1)/mass
    dF[4,1]=(y[1]*F2+2*lambd)/mass
    dF[4,2]=(y[1]*F3)/mass
    dF[4,3]=phiy*y[3]
    dF[4,4]=phiy*y[4]
    dF[4,5]=phiy*y[5]
    dF[5,0]=(y[2]*F1)/mass
    dF[5,1]=(y[2]*F2)/mass
    dF[5,2]=(y[2]*F3+2*lambd)/mass
    dF[5,3]=phiz*y[3]
    dF[5,4]=phiz*y[4]
    dF[5,5]=phiz*y[5]

    return dF

[我正在尝试建立一个使用函数和dF的隐式求解器，以便猜测值，以便为我的摆函数创建更准确的求解器。

我得到的错误如下

追踪（最近通话）：第206行，在test_function（backwards_Euler）中第186行，在test_function y1 = test_function（func，y0，t）中第166行，在向后的Euler中，F = np.asarray（y [i] + dt * function（zold，time [i + 1]）-zold）第13行，以func lambd =（grav.dot（x）+ v.dot（v））/ x.dot（x）

[ValueError：形状（3,6）和（3,6）不对齐：6（dim 1）！= 3（dim 0）] >>

我认为这是因为我正在进行的一种计算正在趋于无穷大，但我无法弄清楚可能导致此问题的原因。

这是我当前的向后欧拉函数：

### Backwards Euler Method

def backwards_Euler(function, y_matrix, time):
    y = np.zeros((np.size(time), np.size(y_matrix)))
    y[0, :] = y_matrix

    dt = time[1] - time[0]

    for i in range(len(time-1)):
        err = 1
        zold = y[i] + dt*function(y[i],time[i])  ## guess with forward euler
        I = 0

        while err > 10**(-10) and I < 5:
            F = y[i] + dt * function(zold, time[i+1])-zold  ## Here is where my error occurs
            dF = dt*dF_matrix(y[i+1])-1
            znew = zold - F/dF
            zold = znew
            I+=1
        y[i+1]=znew
    return y
这是我目前如何测试功能的方法，尽管目前我没有任何输出

# initial condition
y0 = np.array([0.0, 1.0, 0.0, 0.8, 0.0, 1.2])


def test_function(test_function):
    print(test_function.__name__ + "...")
    nt = 2500
    time_start = process_time()
    # time points
    t = np.linspace(0, 25, nt)
    # solve ODE
    y1 = test_function(func, y0, t)
    time_elapsed = (process_time() - time_start)
    print('elapsed time', time_elapsed)
    # compute residual:
    r = y1[:, 0] ** 2 + y1[:, 1] ** 2 + y1[:, 2] ** 2 - 1
    rmax1 = np.max(np.abs(r))
    print('error', rmax1)

    fig = plt.figure()
    ax = plt.axes(projection='3d')
    ax.plot3D(y1[:, 0], y1[:, 1], y1[:, 2], 'gray')

    plt.show()


test_function(backwards_Euler)

我正在尝试为摆F设置一个隐式求解器，而dF的定义为：def func（y，t）：###简化了消除摩擦的函数，因为它消除了x，v = y [：3 ]，y [3：6] ...

Answer 1

不，错误是某些矩阵运算中的矩阵形状不匹配。直接的原因是，当x.dot(x)是一个简单数组但不适用于二维数组时，x起作用。然后必须严格遵守矩阵乘法的规则。