Python sympy.dsolve函数代码示例

def main():

    print 'Initial metric:'
    pprint(gdd)
    print '-'*40
    print 'Christoffel symbols:'
    for i in [0,1,2,3]:
    	for k in [0,1,2,3]:
    		for l in [0,1,2,3]:
    			if Gamma.udd(i,k,l) != 0 :
		    		pprint_Gamma_udd(i,k,l)
    print'-'*40
    print'Ricci tensor:'
    for i in [0,1,2,3]:
    	for j in [0,1,2,3]:
    		if Rmn.dd(i,j) !=0:
    			pprint_Rmn_dd(i,j)
    
    print '-'*40
    #Solving EFE for A and B
    s = ( Rmn.dd(1,1)/ A(r) ) + ( Rmn.dd(0,0)/ B(r) )
    pprint (s)
    t = dsolve(s, A(r))
    pprint(t)
    metric = gdd.subs(A(r), t)
    print "metric:"
    pprint(metric)
    r22 = Rmn.dd(3,3).subs( A(r), 1/B(r))
    h = dsolve( r22, B(r) )
    pprint(h)

def test_ode10():
    f = Function("f")
    #type:2nd order, constant coefficients (two real different roots)
    eq = Eq(f(x).diff(x,x) - 3*diff(f(x),x) + 2*f(x), 0)
    assert dsolve(eq, f(x)) in [
        Symbol("C1")*exp(2*x) + Symbol("C2")*exp(x),
        Symbol("C1")*exp(x) + Symbol("C2")*exp(2*x),
    ]
    checksol(eq, f(x), dsolve(eq, f(x)))

def edo_main():
    # const = [0]*6 #a0, a1, a2, a3, a4, a5   ; nessa ordem
    print "Valores", const

    ##Criando uma variavel para guardar valores digitados anteriormente,este só vai calcular as coisas novamente se
    # e somente se parâmetros de entrada diferentes forem colocados
    constAnterior = const

    # a5, a4, a3, a2, a1, a0, xT = input_coefs()
    ###Dif equation solver
    ##Sets if it is of homogenous or inhomogenous type and type of resolution method
    if Respostas[7] == 0:
        solvedEq = dsolve(sympify(Respostas[8]), y(t), hint="nth_linear_constant_coeff_homogeneous")
        # elif (a3 != 0) or (a4 != 0) or (a5 != 0):
        # solvedEq = dsolve(sympify(eq),y(t),hint='nth_linear_constant_coeff_variation_of_parameters')

    else:
        solvedEq = dsolve(sympify(Respostas[8]), y(t), hint="nth_linear_constant_coeff_undetermined_coefficients")

    ##Transformação do tipo sympy_unity para o sympy_mul (mais operações permitidas)
    sepEq = solvedEq._args[1]
    sepEq = sepEq.evalf(prec)

    if const[5] != 0:
        RespPart = sepEq.subs([(C1, 0), (C2, 0), (C3, 0), (C4, 0), (C5, 0)])
    elif const[4] != 0:
        RespPart = sepEq.subs([(C1, 0), (C2, 0), (C3, 0), (C4, 0)])
    elif const[3] != 0:
        RespPart = sepEq.subs([(C1, 0), (C2, 0), (C3, 0)])
    elif const[2] != 0:
        RespPart = sepEq.subs([(C1, 0), (C2, 0)])
    elif const[1] != 0:
        RespPart = sepEq.subs(C1, 0)

    ##Resposta transitória alocada em  RespTran, natural em RespNat

    Respostas[3] = RespPart.evalf(prec)
    formaNatural = sepEq.subs(RespPart, 0)
    Respostas[1] = formaNatural  # Adicionando Forma natural de resposta na lista de respostas
    ## fN é a mesma coisa, mas usado por um bug bizarro do Sympy que exige uma variável sem alocações prévias quando diferenciando
    ##isso é válido no método conds_iniciais_aplicadas
    fN = formaNatural
    rP = RespPart.evalf(prec)
    raizes()
    conds_iniciais_aplicadas(fN, rP)

    respForc = Respostas[4] + Respostas[3]  # Yf = Yt + Yp
    Respostas[5] = respForc.evalf(prec)  # Adiciona Resposta Forcada a lista de respostas

    respComp = Respostas[2]  # Resposta completa p/ eqs. homogeneas
    if Respostas[7] != 0:  # Eqs. nao homogeneas
        respComp = Respostas[2] + Respostas[5]  # Respsota completa p/ eqs. nao-homogeneas

    Respostas[6] = respComp.evalf(prec)  # Adiciona Resposta Completa a lista de respostas

    print "Respostas em ODE_main", Respostas

def main():
    x = Symbol("x")
    f = Function("f")

    eq = Eq(f(x).diff(x), f(x))
    print "Solution for ", eq, " : ", dsolve(eq, f(x))

    eq = Eq(f(x).diff(x, 2), -f(x))
    print "Solution for ", eq, " : ", dsolve(eq, f(x))

    eq = Eq(x ** 2 * f(x).diff(x), -3 * x * f(x) + sin(x) / x)
    print "Solution for ", eq, " : ", dsolve(eq, f(x))

def symbolic_solve():
    """
    Solves the problem symbolically using sympy
    """
    f = Function('f')
    sol = dsolve(2 * Derivative(f(x), x, x) - 1, f(x))
    print sol

def separable_equation(g, h, hf = None):
    """
    dy/dx = g(x)*h(y)
    dy/h(y) = g(x)*dx
    """
    dy, dx = symbols('dy, dx')
    print '\nODE to solve:'
    pprint(Eq(dy/dx, g*h))
    pprint(Eq(dy/h, g*dx))
    print '\nintegrate both sides:'
    LHS, RHS = symbols('LHS, RHS')
    pprint(Eq(LHS, Integral(1/h,y)))
    H = integrate(1/h,y)
    pprint(Eq(LHS, H))
    pprint(Eq(RHS,Integral(g,x)))
    G = integrate(g, x)
    pprint(Eq(RHS, G))
    C = symbols('C')
    print '\nsolving LHS = RHS + C...'
    eq = Eq(H,G+C)
    pprint(eq)
    pprint(solve(eq, y))
    if hf:
        print '\nsolving ODE directly ...'
        pprint(dsolve(f(x).diff(x)-g*hf, f(x)))

    def __call__(self, equations, variables=None):
        if variables is None:
            variables = {}

        if equations.is_stochastic:
            raise ValueError('Cannot solve stochastic equations with this state updater')

        diff_eqs = equations.substituted_expressions

        t = Symbol('t', real=True, positive=True)
        dt = Symbol('dt', real=True, positive=True)
        t0 = Symbol('t0', real=True, positive=True)
        f0 = Symbol('f0', real=True)
        # TODO: Shortcut for simple linear equations? Is all this effort really
        #       worth it?

        code = []
        for name, expression in diff_eqs:
            rhs = expression.sympy_expr
            non_constant = _non_constant_symbols(rhs.atoms(),
                                                 variables) - set([name])
            if len(non_constant):
                raise ValueError(('Equation for %s referred to non-constant '
                                  'variables %s') % (name, str(non_constant)))
            # We have to be careful and use the real=True assumption as well,
            # otherwise sympy doesn't consider the symbol a match to the content
            # of the equation
            var = Symbol(name, real=True)
            f = sp.Function(name)
            rhs = rhs.subs(var, f(t))
            derivative = sp.Derivative(f(t), t)
            diff_eq = sp.Eq(derivative, rhs)
            general_solution = sp.dsolve(diff_eq, f(t))
            # Check whether this is an explicit solution
            if not getattr(general_solution, 'lhs', None) == f(t):
                raise ValueError('Cannot explicitly solve: ' + str(diff_eq))
            # Solve for C1 (assuming "var" as the initial value and "t0" as time)
            if Symbol('C1') in general_solution:
                if Symbol('C2') in general_solution:
                    raise ValueError('Too many constants in solution: %s' % str(general_solution))
                constant_solution = sp.solve(general_solution, Symbol('C1'))
                if len(constant_solution) != 1:
                    raise ValueError(("Couldn't solve for the constant "
                                      "C1 in : %s ") % str(general_solution))
                constant = constant_solution[0].subs(t, t0).subs(f(t0), var)
                solution = general_solution.rhs.subs('C1', constant)
            else:
                solution = general_solution.rhs.subs(t, t0).subs(f(t0), var)
            # Evaluate the expression for one timestep
            solution = solution.subs(t, t + dt).subs(t0, t)
            # only try symplifying it -- it sometimes raises an error
            try:
                solution = solution.simplify()
            except ValueError:
                pass

            code.append(name + ' = ' + sympy_to_str(solution))

        return '\n'.join(code)

def non_homo_linear(rhs, *cds):
    char_func, eq, rs = homo_linear(*cds)
    print('\nnon-homogeneous ODE:')
    pprint(Eq(eq, rhs))
    eq -= rhs
    print('\nsolving non-homogeneous linear equation...\nresult:')
    rs = dsolve(eq)
    pprint(rs)
    return char_func, eq, rs

def test_make_ode_01():
    ode, params = _make_ode_01()
    t, y, y0, k = params
    result = dsolve(ode, y[1](t))
    eq_assumption = sympy.Q.is_true(Eq(k[1], k[0]))
    refined = result.refine(~eq_assumption)
    ignore = k + y0 + (t,)
    int_const = [fs for fs in refined.free_symbols if fs not in ignore][0]
    ref = int_const*exp(-k[1]*t) - exp(-k[0]*t)*k[0]*y0[0]/(k[0] - k[1])
    assert (refined.rhs - ref).simplify() == 0

    def solve_model(self, ics, d_ics):
        model_soln = [dsolve(eq) for eq in self.model]

        def solve_constants(eq, ics, d_ics):
            udiff = Eq(d_ics[0][1], eq.rhs.diff(t))
            system = [eq.subs(ics), udiff.subs(t, 0)]
            consts = solve(system, [C1, C2])
            return eq.subs(consts)

        model_soln_f = [solve_constants(eq[1], ics[eq[0]], d_ics[eq[0]])
                        for eq in enumerate(model_soln[:len(ics)])]

        self.model_response = model_soln_f
        self.x = lambdify(t, model_soln_f[0].rhs, 'numpy')
        self.y = lambdify(t, model_soln_f[1].rhs, 'numpy')

def codegen(expr, lang, indent='    ', ics=None):
    if lang == 'C':
        code = ccode
    elif lang == 'Fortran':
        code = fcode
    else:
        raise ValueError("Lang must be 'C' or 'Fortran'")

    try:
        sol = dsolve(expr, ics=ics)
    except ValueError:
        # Not an ODE
        return code(expr)

    return ccode(sol.rhs, assign_to=sol.lhs.func.__name__)

def solveSystem(L,coordinates,initialConditions):
    """
    This is currently not in use; I'm not sure it ever will be used.    
    
    L is a sympy expression for the Lagrangian.
    coordinates is a list of tuples, like so:
    [(x,xdot),(y,ydot)]
    initialConditions are a list of the above form.
    """    
    (eulerLagrange,coordinates_t,t)=calcEL(L, coordinates)
    for (i,(EL_i,coordinates_t_i)) in\
            enumerate(zip(eulerLagrange,coordinates_t)):
        eqn=sp.dsolve(EL_i,coordinates_t_i(t))
        freeVars=filter(lambda x:x!=t,eqn.atoms())
        newFreeVars=[sp.Symbol(str(freeVar)+"_%i"%i)
            for freeVar in freeVars]

def homo_linear(*cds):
    char_func = 0
    eq = 0
    for cd in cds:
        try:
            c, d = cd
        except TypeError:
            c, d = 1, cd
        char_func += c * y**d
        eq += c * f(x).diff(x, d)
    print('\nhomogeneous ODE:')
    pprint(Eq(eq, 0))
    print('\nhomogeneous characteristic function:')
    pprint(char_func)
    print('\nsolving characteristic function...\nresult:')
    pprint(solve(char_func))
    print('\nsolving homogeneous linear equation...\nresult:')
    rs = dsolve(eq)
    pprint(rs)
    return char_func, eq, rs

def main():

    print("Initial metric:")
    pprint(gdd)
    print("-"*40)
    print("Christoffel symbols:")
    pprint_Gamma_udd(0, 1, 0)
    pprint_Gamma_udd(0, 0, 1)
    print()
    pprint_Gamma_udd(1, 0, 0)
    pprint_Gamma_udd(1, 1, 1)
    pprint_Gamma_udd(1, 2, 2)
    pprint_Gamma_udd(1, 3, 3)
    print()
    pprint_Gamma_udd(2, 2, 1)
    pprint_Gamma_udd(2, 1, 2)
    pprint_Gamma_udd(2, 3, 3)
    print()
    pprint_Gamma_udd(3, 2, 3)
    pprint_Gamma_udd(3, 3, 2)
    pprint_Gamma_udd(3, 1, 3)
    pprint_Gamma_udd(3, 3, 1)
    print("-"*40)
    print("Ricci tensor:")
    pprint_Rmn_dd(0, 0)
    e = Rmn.dd(1, 1)
    pprint_Rmn_dd(1, 1)
    pprint_Rmn_dd(2, 2)
    pprint_Rmn_dd(3, 3)
    # print()
    # print "scalar curvature:"
    # print curvature(Rmn)
    print("-"*40)
    print("Solve Einstein's equations:")
    e = e.subs(nu(r), -lam(r)).doit()
    l = dsolve(e, lam(r))
    pprint(l)
    lamsol = solve(l, lam(r))[0]
    metric = gdd.subs(lam(r), lamsol).subs(nu(r), -lamsol)  # .combine()
    print("metric:")
    pprint(metric)

 def process(self):
     """
     Procesamos la ecuacion y obtenemos los resultados
     """
     
     # Guardamos expresion
     expr = ""
     
     # Obtenemos los symbolos
     x = sp.Symbol("x")
     y = sp.Function("y")
     
     # Obtenemos la expresion
     ec = sp.sympify(self.ec)
     
     # Valor inicial (PVI)
     pvi = { y(self.x): self.y }
     
     # preparamos la EDO
     edo = sp.Eq(y(x).diff(x), ec)
     expr += "EDO:\n\t"+str(edo)
     
     # Despejamos Y
     res = sp.dsolve(y(x).diff(x) - ec)
     # Obtenemos y(x) = f(x)
     expr += "\nEDO resuelta:\n\t"+str(res)
     
     # reemplazamos PVI
     c_eq = sp.Eq(res.lhs.subs(x, 0).subs(pvi), res.rhs.subs(x, 0))
     expr += "\nRemplazamos PVI:\n\t"+str(c_eq)
     
     # Obtenemos el valor de la constante
     c = sp.solve(c_eq)
     expr += "\nValor de C1:\n\t"+str(c[0])
     
     # almacenamos los valores de interes
     self.c = c[0] # valor de C1
     self.res = res # ecuacion resuelta
     
     # retornamos el resultado
     return expr