Python Matrix.applyfunc方法代码示例

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import applyfunc [as 别名]
def test_applyfunc_matrix():
    double = Lambda(x, x**2)

    expr = ElementwiseApplyFunction(double, Xd)
    assert isinstance(expr, ElementwiseApplyFunction)
    assert expr.doit() == Xd.applyfunc(lambda x: x**2)
    assert expr.shape == (3, 3)

    expr = ElementwiseApplyFunction(double, X)
    assert isinstance(expr, ElementwiseApplyFunction)
    assert isinstance(expr.doit(), ElementwiseApplyFunction)
    assert expr == X.applyfunc(double)

    expr = ElementwiseApplyFunction(exp, X*Y)
    assert expr.expr == X*Y
    assert expr.function == exp
    assert expr == (X*Y).applyfunc(exp)

    assert isinstance(X*expr, MatMul)
    assert (X*expr).shape == (3, 3)
    Z = MatrixSymbol("Z", 2, 3)
    assert (Z*expr).shape == (2, 3)

    expr = ElementwiseApplyFunction(exp, Z.T)*ElementwiseApplyFunction(exp, Z)
    assert expr.shape == (3, 3)
    expr = ElementwiseApplyFunction(exp, Z)*ElementwiseApplyFunction(exp, Z.T)
    assert expr.shape == (2, 2)

    raises(ShapeError, lambda: ElementwiseApplyFunction(exp, Z)*ElementwiseApplyFunction(exp, Z))

    M = Matrix([[x, y], [z, t]])
    expr = ElementwiseApplyFunction(sin, M)
    assert isinstance(expr, ElementwiseApplyFunction)
    assert expr.function == sin
    assert expr.expr == M
    assert expr.doit() == M.applyfunc(sin)
    assert expr.doit() == Matrix([[sin(x), sin(y)], [sin(z), sin(t)]])

    expr = ElementwiseApplyFunction(double, Xk)
    assert expr.doit() == expr
    assert expr.subs(k, 2).shape == (2, 2)
    assert (expr*expr).shape == (k, k)
    M = MatrixSymbol("M", k, t)
    expr2 = M.T*expr*M
    assert isinstance(expr2, MatMul)
    assert expr2.args[1] == expr
    assert expr2.shape == (t, t)
    expr3 = expr*M
    assert expr3.shape == (k, t)

    raises(ShapeError, lambda: M*expr)

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import applyfunc [as 别名]
def _jac(robo, symo, n, i, j, chain=None, forced=False, trig_subs=False):
    """
    Computes jacobian of frame n (with origin On in Oj) projected to frame i
    """
    #  symo.write_geom_param(robo, 'Jacobian')
    # TODO: Check projection frames, rewrite DGM call for higher efficiency
    J_col_list = []
    if chain is None:
        chain = robo.chain(n)
        chain.reverse()
    # chain_ext = chain + [robo.ant[min(chain)]]
    # if not i in chain_ext:
    #     i = min(chain_ext)
    # if not j in chain_ext:
    #     j = max(chain_ext)
    kTj_dict = dgm(robo, symo, chain[0], j, key='left', trig_subs=trig_subs)
    kTj_tmp = dgm(robo, symo, chain[-1], j, key='left', trig_subs=trig_subs)
    kTj_dict.update(kTj_tmp)
    iTk_dict = dgm(robo, symo, i, chain[0], key='right', trig_subs=trig_subs)
    iTk_tmp = dgm(robo, symo, i, chain[-1], key='right', trig_subs=trig_subs)
    iTk_dict.update(iTk_tmp)
    for k in chain:
        kTj = kTj_dict[k, j]
        iTk = iTk_dict[i, k]
        isk, ink, iak = Transform.sna(iTk)
        sigm = robo.sigma[k]
        if sigm == 1:
            dvdq = iak
            J_col = dvdq.col_join(Matrix([0, 0, 0]))
        elif sigm == 0:
            dvdq = kTj[0, 3]*ink-kTj[1, 3]*isk
            J_col = dvdq.col_join(iak)
        else:
            J_col = Matrix([0, 0, 0, 0, 0, 0])
        J_col_list.append(J_col.T)
    Jac = Matrix(J_col_list).T
    Jac = Jac.applyfunc(symo.simp)
    iRj = Transform.R(iTk_dict[i, j])
    jTn = dgm(robo, symo, j, n, fast_form=False, trig_subs=trig_subs)
    jPn = Transform.P(jTn)
    L = -tools.skew(iRj*jPn)
    L = L.applyfunc(trigsimp)
    if forced:
        symo.mat_replace(Jac, 'J', '', forced)
        L = symo.mat_replace(L, 'L', '', forced)
    return Jac, L

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import applyfunc [as 别名]

#.........这里部分代码省略.........
        Vector.basis = tuple(basis)

        # define metric tensor

        default_metric = []
        for bv1 in Vector.basis:
            row = []
            for bv2 in Vector.basis:
                row.append(Vector.basic_dot(bv1, bv2))
            default_metric.append(row)
        Vector.metric = Matrix(default_metric)
        if metric is not None:
            if metric[0] == '[' and metric[-1] == ']':
                Vector.is_orthogonal = True
                metric_str_lst = metric[1:-1].split(',')
                Vector.metric = []
                for g_ii in metric_str_lst:
                    Vector.metric.append(S(g_ii))
                Vector.metric = Matrix(Vector.metric)
            else:
                metric_str_lst = flatten(str_array(metric))
                for index in range(len(metric_str_lst)):
                    if metric_str_lst[index] != '#':
                        Vector.metric[index] = S(metric_str_lst[index])

        Vector.metric_dict = {}  # Used to calculate dot product
        N = range(len(Vector.basis))
        if Vector.is_orthogonal:
            for ii in N:
                    Vector.metric_dict[Vector.basis[ii].obj] = Vector.metric[ii]
        else:
            for irow in N:
                for icol in N:
                    Vector.metric_dict[(Vector.basis[irow].obj, Vector.basis[icol].obj)] = Vector.metric[irow, icol]

        # calculate tangent vectors and metric for curvilinear basis

        if curv != (None, None):
            X = S.Zero
            for (coef, base) in zip(curv[0], Vector.basis):
                X += coef * base.obj
            Vector.tangents = []
            for (coord, norm) in zip(Vector.coords, curv[1]):
                tau = diff(X, coord)
                tau = trigsimp(tau)
                tau /= norm
                tau = expand(tau)
                Vtau = Vector()
                Vtau.obj = tau
                Vector.tangents.append(Vtau)
            metric = []
            for tv1 in Vector.tangents:
                row = []
                for tv2 in Vector.tangents:
                    row.append(tv1 * tv2)
                metric.append(row)
            metric = Matrix(metric)
            metric = metric.applyfunc(TrigSimp)
            Vector.metric_dict = {}
            if metric.is_diagonal:
                Vector.is_orthogonal = True
                tmp_metric = []
                for ii in N:
                    tmp_metric.append(metric[ii, ii])
                    Vector.metric_dict[Vector.basis[ii].obj] = metric[ii, ii]
                Vector.metric = Matrix(tmp_metric)
            else:
                Vector.is_orthogonal = False
                Vector.metric = metric
                for irow in N:
                    for icol in N:
                        Vector.metric_dict[(Vector.basis[irow].obj, Vector.basis[icol].obj)] = Vector.metric[irow, icol]
            Vector.norm = curv[1]

            if debug:
                oprint('Tangent Vectors', Vector.tangents,
                       'Metric', Vector.metric,
                       'Metric Dictionary', Vector.metric_dict,
                       'Normalization', Vector.norm, dict_mode=True)

            # calculate derivatives of tangent vectors

            Vector.dtau_dict = None
            dtau_dict = {}

            for x in Vector.coords:
                for (tau, base) in zip(Vector.tangents, Vector.basis):
                    dtau = tau.diff(x).applyfunc(TrigSimp)
                    result = S.Zero
                    for (t, b) in zip(Vector.tangents, Vector.basis):
                        t_dtau = TrigSimp(t * dtau)
                        result += t_dtau * b.obj
                    dtau_dict[(base.obj, x)] = result

            Vector.dtau_dict = dtau_dict

            if debug:
                oprint('Basis Derivatives', Vector.dtau_dict, dict_mode=True)

        return tuple(Vector.basis)

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import applyfunc [as 别名]
def Jump_H(x,x_i,Xi):
    return 1./2.*(SGN_H(x-Xi)-SGN_H(x_i-Xi))

N_mtx_S = Matrix([1-x,x, (1-x)*Jump_S(x,0.,Xi_jump), x*Jump_S(x,1.,Xi_jump)])
N_mtx_H = Matrix([1-x,x, (1-x)*Jump_H(x,0.,Xi_jump), x*Jump_H(x,1.,Xi_jump)])

print "N_s\n",N_mtx_S
print "N_h\n",N_mtx_H

dx = lambda a: diff(a,x)
p = Plot()
for i,N in enumerate(N_mtx_S):
    p[i]=N,[x,0,1,100]
    
B_mtx_S = N_mtx_S.applyfunc(dx)
B_mtx_H = N_mtx_H.applyfunc(dx)

print 'B_s\n',B_mtx_S
print 'B_h\n',B_mtx_H

Ef, Em = symbols('Ef,Em') 
D_mtx = Matrix([Ef + Em *(1 - DiracDelta(Xi_jump))])

K_r_S = B_mtx_S * D_mtx*B_mtx_S.T
K_r_H = B_mtx_H * D_mtx*B_mtx_H.T

K_r_Sdata = K_r_S.subs(Ef,20).subs(Em,20)
K_r_Hdata = K_r_H.subs(Ef,20).subs(Em,20)

print 'K_r_s\n',K_r_S.shape

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import applyfunc [as 别名]
    return 1./2.*(SGN_H(x-Xi)-SGN_H(x_i-Xi))

N_mtx_S = Matrix(1,4,[(x-1)*(y-1),-x*(y-1), (x-1)*(y-1)*Jump_S(x,0.,Xi_jump), -x*(y-1)*Jump_S(x,1.,Xi_jump)])
N_mtx_H = Matrix(1,4,[(x-1)*(y-1),-x*(y-1), (x-1)*(y-1)*Jump_H(x,0.,Xi_jump), -x*(y-1)*Jump_H(x,1.,Xi_jump)])

print "N_s\n",N_mtx_S
print "N_h\n",N_mtx_H

dx = lambda a: diff(a,x)
dy = lambda a: diff(a,y)

#p = Plot()
#for i,N in enumerate(N_mtx_S):
#    p[i]=N,[x,0,1,100]
    
B_mtx_S_x = N_mtx_S.applyfunc(dx)
B_mtx_H_x = N_mtx_H.applyfunc(dx)

#print 'B_sx\n',B_mtx_S_x
#print 'B_hx\n',B_mtx_H_x

B_mtx_S_xy = N_mtx_S.applyfunc(dy)
B_mtx_H_xy = N_mtx_H.applyfunc(dy)

#print 'B_sxy\n',B_mtx_S_xy
#print 'B_hxy\n',B_mtx_H_xy

B_mtx_S = B_mtx_S_x.col_join(B_mtx_S_xy)
B_mtx_H = B_mtx_H_x.col_join(B_mtx_H_xy)

print 'B_s\n',B_mtx_S