Code for Higher Mathematics I First Order ODEs
Table of Contents
These are the code snippets used in First Order ODEs part of Higher Mathematics I.
Introduction
All code herein are solved using SAGEMath, a python based library focusing on mathematics.
Solve the initial value problem
\[\rp{\cos y \sinh x} + 1 \dd\x sin\dd \y \
# Define variables and functions x = var('x') y = function('y')(x) # Define function func = diff(y,x) - (cos(y)*sinh(x) + 1) / (sin(y) * cosh(x)) == 0 # Solve Equation desolve(func, y, ics=[1,2])
-cos(y(x))*cosh(x) - x == -cos(2)*cosh(1) - 1
x, y = var('x, y') func = -cos(y)*cosh(x) - x == -cos(2)*cosh(1) - 1 plt = implicit_plot(func,(x , 0, 3) ,(y, -1, 3),linewidth=2, color='purple') kwargs = { 'xmin' : -0.5 , 'xmax' : +3 , 'ymin' : -0.5 , 'ymax' : +3 , 'gridlines' : 'minor' , 'frame' : False , 'axes' : True , 'fig_tight': True, 'gridlinesstyle': dict(color="#E5E4E2", linestyle=":", linewidth=0.25), 'transparent' : True, 'figsize' : 6, 'aspect_ratio' : 7/12, } plt.save("images/First-Order-ODEs/ivp-ex-a.pdf",**kwargs)
Example: Ötzi the Iceman
In September 1991 the famous Iceman (Oetzi), a mummy from the Neolithic period of the Stone Age found in the ice of the Oetztal Alps (hence the name “Oetzi”) in Southern Tyrolia near the Austrian–Italian border, caused a scientific sensation.
When did Oetzi approximately live and die if the ratio of carbon 164C to carbon 162C in this mummy is 52.5% of that of a living organism? Physical Information. In the atmosphere and in living organisms, the ratio of radioactive carbon 164C (made radioactive by cosmic rays) to ordinary carbon 162C is constant. When an organism dies, its absorption of 164C by breathing and eating terminates.
Hence one can estimate the age of a fossil by comparing the radioactive carbon ratio in the fossil with that in the atmosphere. To do this, one needs to know the half-life of 164C, which is 5715 years (CRC Handbook of Chemistry and Physics, 83rd ed., Boca Raton: CRC Press, 2002, page 11–52, line 9).
# Define variables and functions k, t, _C = var('k, t, _C') y = function('y')(t) # Define function func = diff(y, t) - k * y == 0 # Solve Equation eq = desolve(func, y, ivar = t) # Use the half-life value of t = 5715 sol = solve(eq == 0.5*_C, k) _inter = sol[0].subs(t = 5715) simplify(solve(e**(_inter.rhs()*t) == 0.525, t)[0].rhs().n())
5312.72499110066
It seems our iceman was around 5300 years ago.
Exercise: The Mixing Problem
Mixing problems occur quite frequently in chemical industry. We explain here how to solve the basic model involving a single tank. The tank in Fig. 11 contains 1000 gal of water in which initially 100 lb of salt is dissolved. Brine runs in at a rate of 10 gal>min, and each gallon contains 5 lb of dissoved salt.
The mixture in the tank is kept uniform by stirring. Brine runs out at 10 gal>min. Find the amount of salt in the tank at any time t.
# Define variables and functions A, B, t = var('A, B, t') y = function('y')(t) # Define function func = diff(y, t) - A + B * y == 0 # Solve Equation eq = desolve(func, y, ics = [0, 100], ivar = t) _inter = eq.subs(A = 50, B = 0.01) res = _inter.plot(t, 0, 500) res.save("mixing_problem.pdf")
Exercise 2
Using Theorem 1 or 2, find an integrating factor and solve the initial value problem.
x = var('x') y = function('y')(x) func = (diff(y,x) - exp(x + y) + y * exp(y)) / (1 - x * exp(y)) == 0 desolve(func, y, ics=[0,-1])
# Define the variables A, B, K, w, t = var('A, B, K, w, t') y = function('y')(t) # Define the function func = diff(y, t) + K * y - A + B * cos (w * t) == 0 # Solve the linear equation hormone = desolve(func, y, ivar = t, ics=[0,0])
# Enter the substitution variables ploet = test.subs(A=1, B=1, K=0.05, w = pi/12) plt = ploet.plot(t,0,200) kwargs = { 'xmin' : 0 , 'xmax' : +200 , 'ymin' : 0 , 'ymax' : +25 , 'gridlines' : 'minor' , 'frame' : False , 'axes' : True , 'fig_tight': True, 'gridlinesstyle': dict(color="#E5E4E2", linestyle=":", linewidth=0.25), 'transparent' : True, 'figsize' : 6, 'aspect_ratio' : 10/3, } plt.save("images/First-Order-ODEs/hormone-levels.pdf",**kwargs)
57/10*e^(3*x)
A= var('A') B= var('B') t= var('t') y = function('y')(t) func = diff(y, t) == A * y - B * y ** 2 desolve(func, y, ivar = t)
x,y = var('x, y') y = function('y')(x) func = diff(y, x) == (e**(x + y) + y * e ** y) / ( 1 - x * e** y) desolve(func, y, ivar = x, ics = [0 , -1])
((x*y(x) + e^x)*e^y(x) + 1)*e^(-y(x)) == e + 1
Time to plot it:
x, y = var('x, y') func = ((x*y + e^x)*e^y + 1)*e^(-y) == e + 1 p = implicit_plot(func,(x , -5, 2) ,(y, -3, 3)) p.save("test.pdf", gridlines=True)
Example: An Exact ODE
Solve the following equation
\begin{equation} \cos \rp{x + y} \dd x + \rp{3x^2 + 2y + \cos \rp{x + y}} \dd y = 0 \end{equation}x = var('x') y = function('y')(x) func = cos(x + y) + (3*y**2 + 2*y + cos(x + y))* diff(y, x) == 0 desolve(func, y, ivar = x)
y(x)^3 + y(x)^2 + sin(x + y(x)) == _C
where _C is an arbitrary constant.
Example: Simple Calculus
Solve the following ODEs by integration.
\begin{align*} y' &= - \sin \pi x, \\ y' &= e ^{-3x}, \\ y' &= x e^{x^2/2}, \\ y' &= \cosh{4x}. \end{align*}from sage.symbolic.integration.integral import indefinite_integral x = var('x') y = function('y')(x) print("(i) ", indefinite_integral(-sin(pi*x), x)) print("(ii) ", indefinite_integral(e**(-3*x), x)) print("(iii) ", indefinite_integral(x*e**(x**2/2), x)) print("(iv) ", indefinite_integral(cosh(4*x), x))
(i) cos(pi*x)/pi (ii) -1/3*e^(-3*x) (iii) e^(1/2*x^2) (iv) 1/4*sinh(4*x)
Example: Subsonic Flight
The efficiency of the engines of subsonic airplanes depends on air pressure and usually is maximum near about 36 000 ft.
Find the air pressure \(y(x)\) at this height without calculation.
Note: The rate of change (y'(x)) is proportional to the pressure, and at 18 000 ft the pressure has decreased to half its value \(y_0\) at sea level.
x, k, _C = var('x, k, _C') y = function('y')(x) func = diff(y, x) - k * y == 0 _inter = desolve(func, y, ivar = x) solve(_inter.subs(x = 18000) - _C * 1/2 == 0, k)
_C*e^(18000*k)
Example: First-Order ODE, General Solution, Initial Value Problem
Solve the initial problem
\[y' + y \tan x = \sin 2x, \qquad y(0) = 1 \]
# Define variables and functions x = var('x') y = function('y')(x) func = diff(y, x) + y * tan(x) == sin(2*x) desolve(func, y, ivar= x, ics = [0 , 1])
-(2*cos(x) - 3)/sec(x)
Example: Electric Circuit
R, L, E, t = var('R, L, E, t') I = function('I')(t) func = diff(I, t) + R / L * I == E / L _inter = desolve(func, I, ivar = t, ics = [0, 0]) eq = simplify(_inter.subs(E=48, R=11, L=0.1)) print(eq)
48/11*(e^(110.0*t) - 1)*e^(-110.0*t)
plt = _inter.subs(E=48, R=11, L=0.1).plot(t,0,0.05) kwargs = { 'xmin' : 0 , 'xmax' : +0.05 , 'ymin' : 0 , 'ymax' : +8 , 'gridlines' : 'minor' , 'frame' : False , 'axes' : True , 'fig_tight': True, 'gridlinesstyle': dict(color="#E5E4E2", linestyle=":", linewidth=0.25), 'transparent' : True, } plt.save("images/First-Order-ODEs/electrical-circuit.pdf",**kwargs)
x = var('x') y = function('y')(x) # Define function func = diff(y, x) - 3 * y == 0 # Solve Equation desolve(func, y, ivar = x, ics=[0, 5.7])
57/10*e^(3*x)
x = var('x') y = function('y')(x) # Define function func = diff(y,x) + y * tan(x) == sin(x) # Solve Equation desolve(func, y, ivar = x)
1/2*(2*_C - log(-sin(x)^2 + 1))/sec(x)
x = var('x') y = function('y')(x) # Define function func = diff(y, x) - cos(x) == 0 # Solve Equation desolve(func, y, ivar = x)
_C + sin(x)
x = var('x') y = function('y')(x) # Define function func = diff(y, x) - 3 * y == 0 # Solve Equation desolve(func, y, ivar = x, ics=[0, 5.7])
57/10*e^(3*x)
x = var('x') y = function('y')(x) # Define function func = 2*x*y*diff(y, x) == y**2 - x**2 # Solve Equation eq = solve(desolve(func, y, ivar = x),y)[0] print(eq)
y(x) == -sqrt(-x^2 - x/_C)
x = var('x') y = function('y')(x) # Define function func = diff(y, x) == - cos(x + y) / (3*y**2 + 2*y + cos(x + y)) # Solve Equation eq = desolve(func, y, ivar = x) print(eq)
y(x)^3 + y(x)^2 + sin(x + y(x)) == _C
A, B, t= var('A, B, t') y = function('y')(t) func = diff(y, t) == A * y - B * y ** 2 desolve(func, y, ivar = )