Solve ode with initial conditions

[jmarcell...@ufpi.edu.br]

I want to solve an ode with initial conditions. For example
using SymPy
f   = SymFunction("f")
x,y = Sym("x,y")

edo1 = SymPy.diff(f(x),x) + x

SymPy.dsolve(edo1)

f(x) == C1 - x^2/2

How to solve C1?

[j verzani]

There may be a better way, but here is one that works:

* extract the RHS of the output: eqn = rhs(out)

* find the symbol for C1 from the output of free_symbols(eqn)

* if (x0,y0) are your initial conditions, then: solve(eqn(x=>x0) - y0, C1)

[jmarcell...@ufpi.edu.br]

Hi j verzani


 SymPy documentation has the "ics", but does not show example: dsolve(eq, func=None, hint='default', simplify=True, ics=None, xi=None, eta=None, x0=0, n=6, **kwargs)
any tips?

[j verzani]

Then this should work where `u` is the Symbolic Function and eqn is just `diff(u(x),x) - 2u(x)`: `dsolve(eqn, ics=Dict(u(x0)=>y0))`, but it doesn't for me :( I think the problem is that initial conditions seem to only work with the power series method. For this problem this call will work with an initial condition: `dsolve(eqn, ics=Dict(u(x0)=>y0), hint="1st_power_series")`. Which isn't too satisfying, as there is an exact answer.

--J

[jmarcell...@ufpi.edu.br]

Hi j
in fact, the solution is not easy. I got your tip and made this way:

C1 = Sym("C1")

# eq1 na forma: substituir x por 0 na parte direita da sol_simb
# e subtrair a expressao pelo valor de y0

eq1 = subs(rhs(sol_simb),x=>0) - 1

SymPy.solve(eq1,C1)

[j verzani]

It can be a bit cleaner to just get the constant if you substitute as the following:

eqn = diff(u(x),x) - 2u(x)   # for example

out = dsolve(eqn)

solve(out(u(x) => y0, x=> x0))

[jmarcell...@ufpi.edu.br]

Hi j
very good,thanks! very practical. For a second order ode, as would be?

[j verzani]

Similar things can be done. Rather than illustrate, I'll point out instead that I just added an interface to `dsolve` that makes this much more straightforward (https://github.com/jverzani/SymPy.jl/blob/master/src/dsolve.jl#L87). I'll tag a new version soon.

[jmarcell...@ufpi.edu.br]

Very good. I hope looking forward to the new version. Idea Code:

using SymPy

y   = SymFunction("y")
x   = Sym("x")

edo2 = SymPy.diff(y(x),x,x) - 4*y(x)-x

sol_simb_ord2 = SymPy.dsolve(edo2)

C1,C2 = Sym("C1,C2")

eq1 = sol_simb_ord2(y(x) => 1, x=> 0)

eq2 = SymPy.diff((sol_simb_ord2),x)(diff(y(x)) => 0, x=> 0)

coef = solve([eq1,eq2])

sol_geral = SymPy.subs(sol_simb_ord2,C1=>coef[C1],C2=>coef[C2])

g(x) = rhs(sol_geral)

g(x)

[j verzani]

Yes, this now can be done with:

y = IVPSolution("y")

x  = Sym("x")

eqn = y''(x) - 4y(x) - x

ivpsolve(eqn, x, (y, 0, 1), (y', 0, 0))

I added a few examples to the tutorial (https://github.com/jverzani/SymPy.jl/blob/master/examples/tutorial.md#initial-value-problems). I'm not entirely happy with the names and could specify the initial conditions differently `y' => (x0, y0)`, say, but otherwise this seems to simplify this task quite a bit.