To solve simple symbolic differential expressions, add the Calculus.jl package to your chosen software.
The Calculus package provides tools for working with the basic calculus operations of differentiation and integration. You can use the Calculus package to produce approximate derivatives by several forms of finite differencing or to produce exact derivative using symbolic differentiation. You can also compute definite integrals by different numerical methods.
[In] differentiate("x^2*y", :x)
differentiate("cos(x) + sin(x) + exp(-x) * cos(x)", :x)
differentiate("cos(x) + sin(y) + exp(-x) * cos(y)", [:x, :y])
[Out] :((2 * 1 * x ^ (2 - 1)) * y + x ^ 2 * 0)
:(1 * -(sin(x)) + 1 * cos(x) + ((-1 * exp(-x)) * cos(x) + exp(-x) * (1 * -(sin(x)))))
2-element Array{Any,1}:
:(1 * -(sin(x)) + ((-1 * exp(-x)) * cos(y) + exp(-x) * 0))
:(1 * cos(y) + (0 * cos(y) + exp(-x) * (1 * -(sin(y)))))
As seen, separable equations are approached by moving the "
y
" terms to one side, the "
x
" terms to the other and integrating. This also applies to autonomous equations then.
There are other families of equation types that have exact solutions, and techniques for solution.
Rather than go over these various families, we demonstrate that SymPy
can solve many of these equations symbolically.
Symbolic functions are defined by SymFunction
(or the macro @symfuns):
using CalculusWithJulia # loads `SymPy`, `Plots`
u = SymFunction("u")
We can write our equation by equating it to 0.
u = SymFunction("u")
@vars x a
eqn = u'(x) - a * u(x) * (1 - u(x))
OUT:
Now we can use the dsolve()
function to separate and solve the function
out = dsolve(eqn)
OUT:
To solve this numerically, we define a problem type by giving it the equation, the initial condition, and the timespan to solve over:
using DifferentialEquations
f(u,p,t) = 1.01*u
u0 = 1/2
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan)
After defining a problem, you solve it using solve
. The solvers can be controlled using the available options are described on the Common Solver Options manual page. For example, we can lower the relative tolerance (in order to get a more correct result,
at the cost of more timesteps) by using the command retol
.
With Julia, you can also explicitly choose the algorithm to use. Many of these algorithms are from
recent research and have been shown to be more efficient than the "standard" algorithms. Tsit5()
will
choose the 5th order Tsitouras method. This is the first algorithm to try in most cases.
For our previously defined differential equation, we can use the follow code to produce an solution.
using DifferentialEquations
f(u,p,t) = 1.01*u
u0 = 1/2
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)
The solution to this differential can be vizualized using the following code
using Plots
plot(sol,linewidth=5,title="Solution to the linear ODE",
xaxis="Time (t)",yaxis="u(t)",label="My Solution") # legend=false
1. Consider an electrical circuit containing a capacitor, resistor, and battery. The charge Q(t) on the capacitor satisfies the equation
where R is the (constant) resistance,Cis the (constant) capacitance, and V is the constant voltage supplied by the battery. IfQ(0) = 0, find Q(t) at any time, and sketch the graph of Q versus t.
Solution
using CalculusWithJulia
Q = SymFunction("Q")
@vars t C V R
eqn = R*Q'(t) + Q(t)/C - V
out = dsolve(eqn)
OUT:
Now since Q(0) = 0, the final solution can be reduced to
Plotting the solution with arbitraty numbers for C, V, and R produces the following graph.
using CalculusWithJulia
@vars t
C = 2.5
V = 8
R = 1
Q = C*V*(1-exp(-t/(R*C)));
plot(Q(t),0 ,10 , legend=false)