[Go to the Notes on Diffy Qs home page]

Notes on Diffy Qs - Using Sage for ODE

Press the Evaluate button below to launch each Sage demonstration. You may have to wait a little before the graph or result appears. Be patient. If you want to do something more complicated, see the Sage documentation. If you want to use Sage more and save your work, you can either try installing Sage on your computer (more difficult) or consider using Sage on the cloud at CoCalc.

Plot a function

Here we plot a function. Most common functions are supported, exponents are done with ^ and don't forget to use * when multiplying. Change the first argument to change the function, and in the second argument you can change the range.

xxxxxxxxxx

plot( (x^2+3*x)*sin(x+0.5)-sinh(x), (x,-5,5))

Several functions can be plotted together. Here we plot sine, cosine, and a polynomial at the same time in different colors.

xxxxxxxxxx

plot1 = plot( sin(x), (x,-5,5), color="blue")

plot2 = plot( cos(x), (x,-5,5), color="red")

plot3 = plot( 0.03*(x-3)*x*(x+3), (x,-5,5), color="green")

(plot1+plot2+plot3).show()

Plot a (numerical) solution of a 1st order ODE

Similar as above, but by default only x is declared as a variable, so we have to now declare y. We draw the solution to $y^{\prime }=f(x,y)$ with an initial condition $y(x_0)=y_0.$ You can edit the function and initial conditions by editing the code. The initial condition is changed via the "ics=" line, here you put in "ics=[$x_0$,$y_0$]", we entered "ics=[1,2]" below meaning $y(1)=2.$ In the "end_points=" you put in where x should start and where it should end, by default $-5\le x\le 5,$ which we entered as "end_points=[-5,5]". If the plot is too jagged, try reducing the step size (see section 1.7, here the method used is Runge-Kutta).

xxxxxxxxxx

var("x y")

f(x,y) = sin(y)+cos(x^2)

sol = desolve_rk4(f(x,y),

                  y, ivar=x,

                  ics=[1,2],

                  end_points=[-5,5],

                  step=0.05, output='plot')

sol.show()

Symbolically solve an ODE

Here, we solve a first order ODE to find a general solution. Here $y^{\prime }$ is written as "diff(y,x)" and the equality must be written as a double equals sign. What we entered in the sample code below is the equation $y^{\prime }=3y+e^x.$ As Sage tries find a solution symbolically here, it will not always succeed. Also note that Sage often gives you an implicit solution (see section 1.3) and leaves figuring out what to do with it up to you.

xxxxxxxxxx

y = function('y')(x)

desolve(diff(y,x) == 3*y + e^x, y).show()

How about an initial condition $y(x_0)=y_0.$ As above, we add "ics=[$x_0$,$y_0$]". In the sample code we used $y(1)=2.$

xxxxxxxxxx

y = function('y')(x)

desolve(diff(y,x) == 3*y + e^x, y, ics=[1,2]).show()

Let's solve second order ODE. Here the second derivative is written as "diff(y,x,2)"

xxxxxxxxxx

y = function('y')(x)

desolve(diff(y,x,2) + diff(y,x) + y == e^x, y).show()

And how about initial conditions for a second order equation. Here we give a triple of numbers for the "ics", that is, $x_0,y(x_0),y^{\prime }(x_0).$ In the sample code, we entered "ics=[0,1,pi]" meaning that $y(0)=1$ and $y^{\prime }(0)=\pi .$

xxxxxxxxxx

y = function('y')(x)

desolve(diff(y,x,2) + diff(y,x) + y == e^x, y, ics=[0,1,pi]).show()

Integrals

Since we're here, we might as well show how to integrate. Let us integrate $\frac{1}{x^2-1}$ first to find the antiderivative and then evaluate the integral from $2$ to $\pi .$ This computation also illustrates that you ought to be careful with computers. The logarithm is not given with absolute value bars as you'd expect, also the logarithm is the natural log but it is output as $\log .$ We again use the "show" function to make things pretty.

xxxxxxxxxx

show(integrate(1/(x^2-1),x))

xxxxxxxxxx

show(integrate(1/(x^2-1),x,2,pi))