UP | HOME

Code for Higher Mathematics I Second Order ODEs

Table of Contents

These are the code snippets used in Second Order ODEs part of Higher Mathematics I.

Introduction

These are the solutions generated by SageMath for the lecture.

Examples

Initial Value Problem

Solve the initial value problem:

\[y\prime\prime + y = 0, \qquad y(0) = 3, \qquad y\prime(0) = -0.5.\] 

x = var('x')
y = function('y')(x)

func_a = diff(y, x, 2) + y == 0

result = desolve(func_a, y, ics = [0,3,-0.5])

print(result)

plt = result.plot(x, 0, 10)

3*cos(x) - 1/2*sin(x)

Reduction of Order if a Solution is Known

Find a basis of solutions of the ODE:

\[ (x^2 - x)y\prime\prime - xy\prime + y = 0 \] 

x, _K1, _K2 = var('x _K1 _K2')
y = function('y')(x)

func_a = diff(y, x, 2) * (x ** 2 - x) - x * diff(y, x, 1) + y  == 0

result = desolve(func_a, y, ivar = x, contrib_ode=True)[0].rhs().subs(_K1 = 1, _K2 = 1)

print(result)

plt = result.plot(x, 0, 10)

x*log(x) + x + 1

BVP: Electric Potential Field Between Two Concentric Spheres

Find the electrostatic potential \(V = v(r)\) between two concentric spheres of radii \(r_1\) = 5 cm and \(r_2\) = 10 cm kept at potentials \(v_1\) = 110 V and \(v_2\) = 0, respectively.

\(V = v(r)\) is a solution of Euler–Cauchy equation \(rv'' + 2v' =\) 0. 

r = var('r')
v = function('v')(r)

func_a = r*diff(v, r, 2) + 2 * diff(v, r, 1) == 0

result = desolve(func_a, v, ics = [5, 110, 10, 0])

print(result)

plt = result.plot(r, 0, 10)

1100/r - 110

IVP: In Case of Distinct Real Roots

Solve the initial value problem

\[y\prime\prime + y\prime - 2y = 0, \quad y(0) = 4, \quad y\prime(0) = -5 \] 

x = var('x')
y = function('y')(x)

func = diff(y, x, 2) + diff(y, x, 1) - 2*y == 0

result = desolve(func, y, ics = [0, 4, -5])

print(result)

plt = result.plot(x, 0, 4)

3*e^(-2*x) + e^x

IVP: In Case of a Double Root

Solve the initial value problem

\[y\prime\prime + y\prime + 0.25y = 0, \quad y(0) = 3, \quad y\prime(0) = -3.5 \] 

x = var('x')
y = function('y')(x)

func = diff(y, x, 2) + diff(y, x, 1) + 0.25*y == 0

result = desolve(func, y, ics = [0, 3, -3.5])

print(result)

plt = result.plot(x, 0, 15)

-(2*x - 3)*e^(-1/2*x)

IVP: In Case of Complex Roots

Solve the initial value problem

\[y\prime\prime + 0.4y\prime + 9.04y = 0, \quad y(0) = 0, \quad y\prime(0) = 3 \] 

x = var('x')
y = function('y')(x)

func = diff(y, x, 2) + 0.4*diff(y, x, 1) + 9.04*y == 0

result = desolve(func, y, ics = [0, 0, 3])

plt = result.plot(x, 0, 30)

Three Cases of Damped Motion

How does the motion of the damped system change if we change the damping constant \(c\) from one to another of the following three values, with:

\begin{align} c &= 100, \\ c &= 60, \\ c &= 10, \end{align}

with \(y(0) = 0.16\) and \(y\prime(0) = 0\). 

x = var('x')
y = function('y')(x)

func_a = 10*diff(y, x, 2) + 100 * diff(y, x, 1) + 90 * y == 0
func_b = 10*diff(y, x, 2) + 60 * diff(y, x, 1) + 90 * y == 0
func_c = 10*diff(y, x, 2) + 10 * diff(y, x, 1) + 90 * y == 0

result_a = desolve(func_a, y, ics = [0, 0.16, 0])
result_b = desolve(func_b, y, ics = [0, 0.16, 0])
result_c = desolve(func_c, y, ics = [0, 0.16, 0])

print(result_a)
print(result_b)
print(result_c)

9/50*e^(-x) - 1/50*e^(-9*x)
4/25*(3*x + 1)*e^(-3*x)
4/875*(sqrt(35)*sin(1/2*sqrt(35)*x) + 35*cos(1/2*sqrt(35)*x))*e^(-1/2*x)

Non-Homogeneous ODEs

Application of the Basic Rule (a)

Solve the initial value problem

\[ y\prime\prime + y = 0.001x^2, \quad y(0) = 0, \quad y\prime(0) = 1.5. \] 

ax = [0 ,40 ,-1.5, 3.5]; ratio = 0.5

kwargs = {
    'xmin' : ax[0] ,
    'xmax' : ax[1] ,
    'ymin' : ax[2] ,
    'ymax' : ax[3] ,
    'gridlines' : 'minor' ,
    'frame' : False ,
    'axes' : True ,
    'fig_tight': True,
    'gridlinesstyle': dict(color="#E5E4E2", linestyle=":", linewidth=0.25),
    'transparent' : True,
    'figsize' : 6,
    'aspect_ratio' : abs(ax[1] - ax[0]) / abs(ax[3] - ax[2]) * ratio,
     }

plt.save("images/Second-Order-ODEs/application-of-the-basic-rule-a.pdf",**kwargs)

Application of the Modification Rule (b)

Solve the initial value problem

\[ y\prime\prime + 3y\prime + 2.25y = -10 e^{-1.5x}, \quad y(0) = 1, \quad y\prime(0) = 0. \] 

x = var('x')
y = function('y')(x)

func = diff(y, x, 2) + 3 * diff(y, x, 1) + 2.25 * y ==  -10 * e^(-1.5 * x)

result = desolve(func, y, ics = [0, 1, 0])

plt = result.plot(x, 0, 5)

Application of the Sum Rule (c)

Solve the initial value problem

\[ y\prime\prime + 2y\prime + 0.75y = 2 \cos x - 0.25 \sin x + 0.09 x, \quad y(0) = 1, \quad y\prime(0) = 0. \] 

x = var('x')
y = function('y')(x)

func = diff(y, x, 2) + 2 * diff(y, x, 1) + 0.75 * y == 2 * cos(x) - 0.25 * sin(x) + 0.09 * x

result = desolve(func, y, ics = [0, 2.78, -0.43])

plt = result.plot(x, 0, 20)

Modelling: Forced Oscillations & Resonance

m, c, k, F_0, w, w_0, t, _K1, _K2 = var('m c k F_0 w w_0 t _K1 _K2')
y = function('y')(t)

c = 0
k = m * w_0 ** 2

assume(w_0 > 0)

func = m * diff(y, t, 2) + c * diff(y, t, 1) + k * y == F_0 * cos(w * t)

result = desolve(func, y, ivar = t).subs(_K1 = 1, _K2 = 1)

rho =  1 / (1 - (w / w_0) ** 2).subs(w = 1)

print("rho is      : ", rho)
print("solution is : ", result)

plt_rho = rho.plot(w_0, 0, 5)

rho is      :  -1/(1/w_0^2 - 1)
solution is :  -F_0*cos(t*w)/(m*w^2 - m*w_0^2) + cos(t*w_0) + sin(t*w_0)

Modelling: Electric Circuits

Find the current \(I(t)\) in an RLC-circuit with:

  • \(R = 11 (ohms)\),
  • \(L = 0.1 (henry)\),
  • \(C = 0.01 (farad)\),

which is connected to a source of:

\[V(t) = 110 \sin (2 \pi 60 t) \]

Assume that current and capacitor charge are 0 when \(t = 0\). 

t = var('t')
y = function('y')(t)

L = 0.1
R = 11
C = 0.01
f = 60
E = 110 * sin ( 2 * pi * f * t)

func = L * diff(y, t, 2) + R * diff(y, t, 1) + 1 / C * y == diff(E, t, 1)

result = desolve(func, y, ivar = t, ics = [0, 0 ,0])

print(result)

plt = result.plot(t, 0, 0.1)

-440/3*pi*e^(-10*t)/(144*pi^2 + 1) + 1100/3*pi*e^(-100*t)/(36*pi^2 + 25) + 660*(66*pi^2*sin(120*pi*t) + (5*pi - 72*pi^3)*cos(120*pi*t))/(5184*pi^4 + 3636*pi^2 + 25)

Variation of Parameters

Method of Variation of Parameters

Solve the non-homogeneous ODE:

\[ y\prime\prime + y = \sec x = \frac{1}{\cos x} \] 

x, _K1, _K2 = var('x _K1 _K2')
y = function('y')(x)

func = diff(y, x, 2) +  y == 1 / cos(x)

result = desolve(func, y).subs(_K1 = 1, _K2 = 1)

print(result)

plt = result.plot(x, 0, 1)

cos(x)*log(cos(x)) + x*sin(x) + cos(x) + sin(x)

Author: Daniel McGuiness