Kernel: SageMath 9.6

C08 Solving First Order Differential Equations with CoCalc

Solving Regular Equations

Exercise 1. Solve the equation ax2+bx+c=0ax^2 + bx + c = 0 for the variable xx in terms of the parameters aa, bb, and cc.

Exercise 2. Solve the equation x7+x1=0x^7 + x - 1 = 0 for xx

Warming or Cooling an Object

Exercise 3. Solve Newton's Law of Cooling differential equation when the ambient temperature is oscillating: u(t)=k(u(t)absin(wt))u'(t) = -k(u(t) - a - b\sin(wt)). How should the variables and parameters be interpreted?

In [0]:
var('t k a b w') #Declare independent variable and parameters u = function('u')(t) #Declare dependent variable de = diff(u, t) == -k*(u - a - b*sin(w*t)); show(de) #Define the ODE usol = desolve(de, dvar = u, ivar = t); show(usol) #Find solution usol = usol.expand(); show(usol) #Simplify solution

Exercise 4. Suppse k=0.2k = 0.2, a=10a = 10, b=5b = 5, and w=1/2w = 1/2. Plot solutions on the interval 0t300 \leq t \leq 30 for C=20C = -20, C=0C = 0, and C=10C = 10

Usain Bolt's Olympic Victory

Table 1.1 reports race splits (seconds) every 10 meters for Usain Bolt's 2008 Olympic gold medal final 100 meter race. The following code draws a scatter plot of this data.

In [0]:
tdata = [0.165, 1.85, 2.87, 3.78, 4.65, 5.50, 6.32, 7.14, 7.96, 8.79, 9.68] pdata = range(0, 110, 10) data = list(zip(tdata, pdata)) dplot = list_plot(data, axes = False, frame = True, axes_labels_size = 1, axes_labels = ["Time (seconds)", "Position (meters)"]) dplot.show(figsize = 4)

Exercise 5. Describe and solve the Hill-Keller differential equation v(t)=Pkv(t)v'(t) = P - kv(t) with initial condition v(t0)=0v(t_0) = 0.

Exercise 6. Choose values for the parameters so that the Hill-Keller model fits Usain Bolt's data well. Plot the data and model on a single set of axes.

Fish Harvesting

Exercise 7. Describe and solve the fish harvesting initial value problem u(t)=ru(t)(1u(t)/K)hu(t),u(0)=u0u'(t) = ru(t)(1-u(t)/K) - hu(t), u(0) = u_0.

Exercise 8. the code below stores Atlantic cod population and harvesting data from Table 1.2 in the text and then shows part of the table. Find the mean of the harvesting data. Graph the population data versus years after January 1, 1978, and the model solution using reasonable parameter values.

In [0]:
ydata = list(range(1978,2009)) udata = [72148, 73793, 74082, 92912, 82323, 59073, 59920, 48789, 70638, 67462, 68702, 61191, 49599, 46266, 34877, 28827, 21980, 17463, 18057, 22681, 20196, 25776, 23796, 19240, 16495, 12167, 21104, 18871, 21241, 22962, 21848] hdata = [0.18847, 0.14974, 0.21921, 0.17678, 0.28203, 0.34528, 0.20655, 0.33819, 0.14724, 0.19757, 0.23154, 0.2086, 0.33565, 0.29534, 0.33185, 0.35039, 0.2827, 0.19928, 0.18781, 0.19357, 0.18953, 0.17011, 0.1566, 0.28179, 0.25287, 0.25542, 0.08103, 0.0874, 0.08195, 0.10518, NaN] table(columns = [ydata[0:32:6], udata[0:32:6], hdata[0:32:6]], header_row = ["year", "$u_t$", "$h_t$"])

Optional code to control the number of digits diplayed for the harvesting data. An error is thrown if NaN enters the formating code.

In [0]:
table(columns = [ydata[1:32:6], udata[1:32:6], list(map(lambda x: '{:.5f}'.format(x), hdata[1:32:6]))], header_row = ["year", "$u_t$", "$h_t$"])

Direction Fields and Graphical Solutions

Exercise 9. Create a direction field for u(t)=0.2(u(t)105sin(0.5t))u'(t) = -0.2(u(t) - 10 - 5\sin(0.5t)).

In [0]:
var('t u') f(t, u) = -0.2*(u - 10 - 5*sin(0.5*t)) #rhs of u' = f(t, u) pd = plot_slope_field(f(t, u), (t, 0, 30), (u, 0, 20)) show(pd, figsize = 4)

Exercise 10. Add many solutions to the direction field.

In [0]:
ps = streamline_plot(f(t, u), (t, 0, 30), (u, 0, 20)) show(pd + ps, figsize = 4)

Exercise 11. Add solution curves for the initial conditions u(0)=0u(0) = 0, u(0)=10u(0) = 10, and u(0)=20u(0) = 20.

Answer 1. Notice that there is no user control over the time interval of each solution. The documentation also states that the initial conditions may be adjusted slightly.

In [0]:
ics = [(0, 0), (0, 10), (0, 20)] ps = streamline_plot(f(t, u), (t, 0, 30), (u, 0, 20), start_points = ics) show(pd + ps, figsize = 4)

Answer 2. This is a more accurate and controllable approach but requires more code..

In [0]:
p1 = desolve_rk4(f(t, u), dvar = u, ivar = t, ics = (0, 0), end_points = 30, output = 'plot') p2 = desolve_rk4(f(t, u), dvar = u, ivar = t, ics = (0, 10), end_points = 30, output = 'plot') p3 = desolve_rk4(f(t, u), dvar = u, ivar = t, ics = (0, 20), end_points = 30, output = 'plot') show(pd + p1 + p2 + p3, figsize = 4)

Answer 3. The supplementary materials to our text have a function that makes the syntax a bit easier. To use this easier syntax, you must have the file draw_dirfield.sage in the same directory as this file and must have rk4_method.sage in the a directory numerics at the same level as the directory in which this file is located. The user does not have control over the size of the figure.

In [0]:
#Load in supplied command for direction fields load('draw_dirfield.sage')
In [0]:
draw_dirfield(f, [0, 30, 0, 20], ics)
In [0]: