The Julia Language

Introduction Part I: Plotting Part II: First Order Differential Equations Part III: Second and Higher Order Differential Equations
Downloading Julia Julia Homepage

First Order Differential Equations

Simple Differentiation

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)))))  

Symbolic Differentiatiom

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:

Using Julia to Solve an ODE

Defining a Problem

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)

Solving a Problem

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)

Plotting a Solution

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

Examples from Class

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)

Exact Equations