Personal Stories

Stories about how you have used Maple, MapleSim and Math in your life or work.

Over the last few days, I’ve been creating worksheets on oscillators to support my class’s understanding of these fundamental physics concepts. I wanted to share one of these worksheets that I found particularly useful for illustrating energy exchange and motion dynamics.

A simple pendulum is a classic physics example that exhibits periodic motion. It consists of a mass m (called a bob) attached to a string or rod of length L, which swings back and forth under the influence of gravity. When the bob is displaced from its equilibrium position and released, it swings back and forth under the influence of gravity.
To derive the equation of motion, we can examine the forces acting on the pendulum bob and use Newton’s second law.

Period of a Pendulum:

• 

Frequency (f) "-" the number of cycles the pendulum completes in one second. Measued in hertz ("Hz)."

f = 1/T

• 

Period ("T) -" the time it takes the pendulum to complete one cycle. Measued in seconds (s).

T = 2*Pi*sqrt(L/g)

This period depends only on the length Land gravitational acceleration "g,"meaning it is independent of the amplitude for small oscillations.

What is the period and the frequency of a single pendulum that is 70 cm long on the earth and on the moon?

L := .7; g__earth := 9.8; g__moon := 1.6

.7

 

9.8

 

1.6

(1)

T__earth := 2*Pi*sqrt(L/g__earth); f__earth := L/T__earth

1.679251909

 

.4168522878

(2)

T__moon := 2*Pi*sqrt(L/g__moon); f__moon := L/T__moon

4.155936442

 

.1684337597

(3)

The above image is taken from https://www.researchgate.net/publication/365297210_Scientific_counterfactuals_as_make-believe

1. 

Forces on the Pendulum Bob:

The main forces acting on the bob are:

• 

The gravitational force"`f__g`=mg, "acting vertically downward.

• 

The tension Tauin the string, acting along the string toward the pivot point.

2. 

Components of the Gravitational Force:

Since the pendulum swings in an arc, it’s helpful to resolve the gravitational force into two components:

• 

Radial Component (along the string): This component, "`f__y`=mgcostheta ," is countered by the tension in the string and does not contribute to the pendulum’s motion.

• 

Tangential Component (perpendicular to the string): This component, f__z = -`mgsinθ`(restoring force), acts along the arc of the pendulum’s swing and is responsible for its motion.

3. 

Applying Newton's Second Law
Since the tangential component of the gravitational force causes the pendulum’s motion, we can apply Newton's second law in the tangential direction:``

f__X = ma__x

Substituting for f__x and the tangential acceleration a__xNULL

m*(diff(s(t), t, t)) = -`mgsinθ`

where diff(s(t), t, t) = a__x and a__x = diff(x, t, t)

Now, we want to write everything in terms of θ

s = `Lθ`

we obtain:

diff(theta(t), t, t) = -g*`sinθ`/L

Small-Angle Approximation

For small angles (typically) theta  , the approximation

`≈`(sin(theta), theta)(where theta is in radians) can be used. This simplifies the equation:

diff(theta(t), t, t) = -g*theta/L

This equation is now in the form of a simple harmonic oscillator

diff(theta(t), t, t) = -omega^2*theta

where omega = sqrt(g/L)is the angular frequency of the pendulum.

restart; with(plots); with(DEtools)

L := 1; m := .2; g := 9.8

1

 

.2

 

9.8

(4)

T := 2*Pi*sqrt(L/g)

2.007089924

(5)

omega := sqrt(g/L)

3.130495168

(6)

ODE__1 := diff(theta(t), t, t)+omega^2*theta = 0; IC := theta(0) = A, (D(theta))(0) = 0

diff(diff(theta(t), t), t)+9.799999997*theta(t) = 0

 

theta(0) = A, (D(theta))(0) = 0

(7)

sol := dsolve({IC, ODE__1}, theta(t))

theta(t) = A*cos((1/100000)*97999999970^(1/2)*t)

(8)

plot_1 := subs(A = 0.873e-1, sol); plotsresult := plot([rhs(plot_1)], t = 0 .. 2, color = [red])

 

`θ_t` := rhs(subs(A = 0.873e-1, sol)); v_t := diff(`θ_t`, t)

0.873e-1*cos((1/100000)*97999999970^(1/2)*t)

 

-0.8730000000e-6*97999999970^(1/2)*sin((1/100000)*97999999970^(1/2)*t)

(9)

T := (1/2)*m*L^2*v_t^2; V := m*g*L*(1-cos(`θ_t`)); H := T+V

0.7468864200e-2*sin((1/100000)*97999999970^(1/2)*t)^2

 

1.96-1.96*cos(0.873e-1*cos((1/100000)*97999999970^(1/2)*t))

 

0.7468864200e-2*sin((1/100000)*97999999970^(1/2)*t)^2+1.96-1.96*cos(0.873e-1*cos((1/100000)*97999999970^(1/2)*t))

(10)

energy_plot := plot([eval(T), eval(V), eval(H)], t = 0 .. 5, color = [red, blue, green], legend = ["Kinetic Energy", "Potential Energy", "Total Energy"], title = "Energy Exchange in Simple Pendulum", labels = ["Time (s)", "Energy (Joules)"])

 

directionfield := DEplot([diff(theta(t), t) = v(t), diff(v(t), t) = -omega^2*theta(t)], [theta(t), v(t)], t = 0 .. 20, theta = -20 .. 20, v = -40 .. 40, arrows = medium, title = "Direction Field for Simple Harmonic Oscillator", axes = boxed, color = navy)

sol1 := dsolve({ODE__1, theta(0) = 3, (D(theta))(0) = 0}, theta(t)); sol2 := dsolve({ODE__1, theta(0) = 6.5, (D(theta))(0) = 0}, theta(t)); sol3 := dsolve({ODE__1, theta(0) = -8, (D(theta))(0) = 0}, theta(t)); sol4 := dsolve({ODE__1, theta(0) = 9.7, (D(theta))(0) = 2.5}, theta(t))

theta(t) = 3*cos((1/100000)*97999999970^(1/2)*t)

 

theta(t) = (13/2)*cos((1/100000)*97999999970^(1/2)*t)

 

theta(t) = -8*cos((1/100000)*97999999970^(1/2)*t)

 

theta(t) = (25000/9799999997)*97999999970^(1/2)*sin((1/100000)*97999999970^(1/2)*t)+(97/10)*cos((1/100000)*97999999970^(1/2)*t)

(11)

theta1 := rhs(sol1); theta2 := rhs(sol2); theta3 := rhs(sol3); theta4 := rhs(sol4)

3*cos((1/100000)*97999999970^(1/2)*t)

 

(13/2)*cos((1/100000)*97999999970^(1/2)*t)

 

-8*cos((1/100000)*97999999970^(1/2)*t)

 

(25000/9799999997)*97999999970^(1/2)*sin((1/100000)*97999999970^(1/2)*t)+(97/10)*cos((1/100000)*97999999970^(1/2)*t)

(12)

v1 := diff(theta1, t); v2 := diff(theta2, t); v3 := diff(theta3, t); v4 := diff(theta4, t)

-(3/100000)*97999999970^(1/2)*sin((1/100000)*97999999970^(1/2)*t)

 

-(13/200000)*97999999970^(1/2)*sin((1/100000)*97999999970^(1/2)*t)

 

(1/12500)*97999999970^(1/2)*sin((1/100000)*97999999970^(1/2)*t)

 

(5/2)*cos((1/100000)*97999999970^(1/2)*t)-(97/1000000)*97999999970^(1/2)*sin((1/100000)*97999999970^(1/2)*t)

(13)

phase_plot := plot([[eval(theta1, t = tval), eval(v1, t = tval), tval = 0 .. 20], [eval(theta2, t = tval), eval(v2, t = tval), tval = 0 .. 20], [eval(theta3, t = tval), eval(v3, t = tval), tval = 0 .. 20], [eval(theta4, t = tval), eval(v4, t = tval), tval = 0 .. 20]], style = line, title = "Phase Portrait for Simple Harmonic Oscillator", labels = ["x (Displacement)", "v (Velocity)"], color = ["red", "blue", "green", "orange"], axes = boxed)

display(directionfield, phase_plot)

 

theta_at_1_sec := subs(t = 1, A = 0.873e-1, rhs(sol)); evalf(theta_at_1_sec)

0.873e-1*cos((1/100000)*97999999970^(1/2))

 

-0.8729462437e-1

(14)
 

NULL

Download Simple_Pendulum.mw

Maple Transactions has just published the Autumn 2024 issue at mapletransactions.org

From the header:

This Autumn Issue contains a "Puzzles" section, with some recherché questions, which we hope you will find to be fun to think about.  The Borwein integral (not the Borwein integral of XKCD fame, another one) set out in that section is, so far as we know, open: we "know" the value of the integral because how could the identity be true for thousands of digits but yet not be really true? Even if there is no proof.  But, Jon and Peter Borwein had this wonderful paper on Strange Series and High Precision Fraud showing examples of just that kind of trickery.  So, we don't know.  Maybe you will be the one to prove it! (Or prove it false.)

We also have some historical papers (one by a student, discussing the work of his great grandfather), and another paper describing what I think is a fun use of Maple not only to compute integrals (and to compute them very rapidly) but which actually required us to make an improvement to a well-known tool in asymptotic evaluation of integrals, namely Watson's Lemma, just to explain why Maple is so successful here.

Finally, we have an important paper on rational interpolation, which tells you how to deal well with interpolation points that are not so well distributed.

Enjoy the issue, and keep your contributions coming.

With the 2024 Maple Conference coming up this week, I imagine one of two of you have noticed something missing. We chose not to have a Conference Art Gallery this year, because we have been working to launch new section of MaplePrimes:  The MaplePrimes Art Gallery. This new section of MaplePrimes is designed for showing off your Maple related images, in a gallery format that puts the images up front, like Instagram but for Math.

To get the ball rolling on the gallery, I have populated it with many of the works from past years' Maple Art Galleries, so you can now browse them all in one place, and maybe give "Thumbs Ups" to art pieces that you felt should have won the coveted "People's Choice Award".

Once you are inspired by past entries, we welcome you to submit your new artwork!  Just like the conference galleries, we are looking for all sorts of mathematical art ranging from computer graphics and animations, to photos of needlework, geometrical sculptures, or almost anything!  Art Gallery posts are very similar to regular MaplePrimes posts except that you are asked to supply an Image File and a Caption that will displayed on the main Gallery Page, the post itself can include a longer description, Maple commands, additional images, and whatever else you like.  See for example this Art Gallery post describing Paul DeMarco's sculpture from the 2021 Maple Conference Gallery - it includes an embedded worksheet as well as additional images.

I can't wait to see what new works of art our MaplePrimes community creates!

 

I recently prepared a worksheet to teach vector fundamentals in one of my classes, and I wanted to share it with you all. It's nothing special, but I found Maple really helpful in demonstrating the concepts visually. Below is a breakdown of what the worksheet covers, with some Maple code examples included.

Feel free to take a look and use it if you find it useful! Any feedback or suggestions on how to improve it would be appreciated.

restart

NULL

v := `<|>`(`<,>`(2, 3)); w := `<|>`(`<,>`(4, 1))

Matrix(%id = 36893488152076804092)

(1)

Basic Vector Operations

• 

Addition and Subtraction

 

We  can add and subtract vectors easily if they are of the same dimension.

NULL

u_add := v+w; u_sub := v-w

Matrix(%id = 36893488152076803596)

(2)

NULL

NULL

Typesetting[delayDotProduct]((((Triangle(L)*a*w*o*f*A*d*d)*i*t*i)*o*n*o)*f, Vector(s), true)

"The famous triangle law can be used for the addition of vectors and this method is also called the head-to-tail method,As per this law,two vectors can be added together by placing them together in such away that the first vector's head joins the tail of the second vector. Thus,by joining the first vector's tail to the head of the second vector,we can obtain the resultant vector sum."

NULL

with(plots); display(arrow([0, 0], [2, 3], color = red, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1), arrow([2, 3], [4, 1], color = green, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1), arrow([0, 0], [6, 4], difference = true, color = blue, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1), scaling = constrained, labels = [x, y], title = "
Triangle Law of Addition of Vectors", titlefont = [times, 20, bold])

 

NULL

NULL

Parallelogram Law of Addition of Vectors

"An other law that can be used for the addition of vectors is the parallelogram law of the addition of vectors*Let's take two vectors v and u,as shown below*They form the two adjacent sides of aparallelogram in their magnitude and direction*The sum v+u is represented in magnitude and direction by the diagonal of the parallelogram through thei rcommon point."

NULL

display(arrow([0, 0], [2, 3], color = red, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1), arrow([0, 0], [6, 4], color = green, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1), arrow([0, 0], [4, 1], difference = true, color = blue, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1), scaling = constrained, labels = [x, y], title = "

Parallelogram Law of Addition of Vectors", titlefont = [times, 20, bold])

 

NULL

NULL

NULL

NULL

NULL

NULL

NULL

• 

Scalar Multiplication

We can multiply a vector by a scalar. To multiply a vector by a scalar (a constant), multiply each of its components by the constant.

Suppose we let the letter  λ represent a real number and  u = (x,y) be the vector

v_scaled := 3*v

Matrix(%id = 36893488152152005076)

(3)

NULL

NULL

• 

-Define*the*opposite*vector*v

Error, (in LinearAlgebra:-Multiply) invalid arguments

 

Two vectors are opposite if they have the same magnitude but opposite direction.

v_opposite := -v; vec1 := arrow([0, 0], [2, 3], color = blue, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1); vec2 := arrow([0, 0], [-2, -3], color = red, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1); display([vec1, vec2], scaling = constrained, title = "Original Vector and its Opposite (-v)")

 
• 

Dot Product

Sometimes the dot product is called the scalar product. The dot product is also an example of an inner product and so on occasion you may hear it called an inner product.

Given the two vectors  `#mover(mi("a"),mo("&rarr;"))` = (x__1, y__1) and `#mover(mi("b"),mo("&rarr;"))` = (x__2, y__2)
 the dot product is, `#mover(mi("a"),mo("&rarr;"))`.`#mover(mi("b"),mo("&rarr;"))` = x__1*x__2+y__1*y__2

`#mover(mi("a"),mo("&rarr;"))`.`#mover(mi("b"),mo("&rarr;"))` = x__1*x__2+y__1*y__2

(4)

`#mover(mi("a"),mo("&rarr;"))` := `<|>`(`<,>`(5, -8))

Matrix(%id = 36893488152255219820)

(5)

`#mover(mi("b"),mo("&rarr;"))` := `<|>`(`<,>`(1, 2))

Matrix(%id = 36893488152255215244)

(6)


display(arrow([0, 0], [5, -8], color = red, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1), arrow([0, 0], [1, 2], difference = true, color = blue, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1), scaling = constrained, labels = [x, y])

 

 

 

 

 

 

 

 

 

 

 

  with(LinearAlgebra)

dot_product := DotProduct(`#mover(mi("a"),mo("&rarr;"))`, `#mover(mi("b"),mo("&rarr;"))`)

-11

(7)
• 

Vector Norm (Magnitude)

To find the magnitude (or length) of a vector, use Norm.

norm_a := Norm(`#mover(mi("a"),mo("&rarr;"))`, 2); norm_b := Norm(`#mover(mi("b"),mo("&rarr;"))`, 2)

5^(1/2)

(8)
• 

Calculate the Cosine Between Two Vectors

cos_theta := dot_product/(norm_a*norm_b)

-(11/445)*89^(1/2)*5^(1/2)

(9)

angle_radians := arccos(cos_theta)

Pi-arccos((11/445)*89^(1/2)*5^(1/2))

(10)

angle_degrees := evalf(convert(angle_radians, degrees))

121.4295656*degrees

(11)

NULL

We can determine whether two vectors are parallel by using the scalar multiple method or the determinant (area of the parallelogram formed by the vectors) method.

``

• 

Scalar Multiple Method

a := Vector([5, -8]); b := Vector([10, -16]); k := a[1]/b[1]; is_parallel := a[2]/b[2] = k

1/2 = 1/2

(12)
• 

Determinant Method

determinant := a[1]*b[2]-a[2]*b[1]; result := determinant = 0

0 = 0

(13)

 

 

display(arrow([0, 0], [5, -8], color = red, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1), arrow([0, 0], [10, -16], difference = true, color = blue, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1), scaling = constrained, labels = [x, y], title = "Parallel Vectors")

 

 

 

 

Two vectors a and b are perpendicular (vertical)

NULL

• 

if and only if their dot product is zero

a1 := Vector([1, 2]); b1 := Vector([-2, 1]); dot_product := DotProduct(a1, b1); is_perpendicular := dot_product = 0

0 = 0

(14)

display(arrow([0, 0], [-2, 1], color = red, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1), arrow([0, 0], [1, 2], difference = true, color = blue, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1), scaling = constrained, labels = [x, y], title = "Perpendicular Vectors")

 
• 

Slope Method

slope_a := a1[2]/a1[1]; slope_b := b1[2]/b1[1]; is_perpendicular := slope_a*slope_b = -1

-1 = -1

(15)

NULL

• 

Special case: If one vector is vertical (undefined slope) and the other is horizontal (zero slope), they are perpendicular.

a2 := Vector([0, 3]); b2 := Vector([4, 0]); is_perpendicular := a2[1] = 0 and b2[2] = 0

true

(16)

vec1 := arrow([0, 0], [4, 0], color = blue, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1); vec2 := arrow([0, 0], [0, 3], color = red, shape = double_arrow, width = 0.1e-1, border = false, head_width = .1, head_length = .1); display([vec1, vec2], scaling = constrained, title = "Perpendicular Vectors")

 

NULL

Download vectors.mw

Just finished an exciting lecture for undergraduate students on Euler methods for solving ordinary differential equations (ODEs)! 📝 We explored the Euler method and the Improved Euler method (a.k.a. Heun's method) and discussed how these fundamental techniques work for approximating solutions to ODEs.

To make things practical and hands-on, I demonstrated how to implement both methods using Maple—a fantastic tool for experimenting with numerical methods. Seeing these concepts come to life with visualizations helps understand the pros and cons of each method.

The Euler method is the simplest numerical approach, where we approximate the solution by stepping along the slope given by the ODE. But there's a catch: using a large step size can lead to large errors that accumulate quickly, causing significant deviation from the real solution.

Plot Results:

  • With a larger step size (h=0.5h = 0.5h=0.5), the solution tends to drift away significantly from the actual curve.
  • Reducing the step size improves accuracy, but it also means more steps and more computation.

Next, we discussed the Improved Euler method, which is like Euler's smarter sibling. Instead of blindly following the initial slope, it:

  1. Estimates the slope at the start.
  2. Uses that slope to predict an intermediate value.
  3. Then, it takes another slope at the intermediate point.
  4. Averages these two slopes for a better approximation.

This technique makes a big difference in accuracy and stability, especially with larger step sizes.

Plot Results:

  • Using the Improved Euler method, we found that even with a larger step size, the solution is smoother and closer to the true path.
  • The average slope helps mitigate the inaccuracies that arise from only using the beginning point's derivative, effectively reducing the local error.
  • The standard Euler method can produce solutions that oscillate or diverge, especially for larger step sizes.
  • The Improved Euler method follows the actual trend of the solution more faithfully, even with fewer steps. This makes it a more efficient choice for balancing computational effort and accuracy.
  • The Euler method is great for its simplicity but often requires very small step sizes to be accurate, leading to more computational effort.
  • The Improved Euler method—which we implemented and visualized on Maple—proved to be more reliable and accurate, especially for larger step sizes.

It was amazing to see the students engage with these foundational methods, and implementing them on Maple brought a deeper understanding of numerical analysis and the challenges of solving differential equations.


 


restart

with(plots)

with(DEtools)

ODE1 := diff(y(x), x) = -2*x*y(x)/(x^2+1)

diff(y(x), x) = -2*x*y(x)/(x^2+1)

(1)

NULL

We calculate the general solution to the ODE

NULL

dsolve(ODE1, y(x))

y(x) = c__1/(x^2+1)

(2)

NULL

Now let's solve the problem for the following two inital conditions

y(0) = 2 and y(0) = 1/2

dsolve({ODE1, y(0) = 2}, y(x))

y(x) = 2/(x^2+1)

(3)

dsolve({ODE1, y(0) = 1/2}, y(x))

y(x) = 1/(2*x^2+2)

(4)

Then, we are going to plot the solutions to the IVP's together with the slope field corresponding to the ODE.

 

NULLdfieldplot(ODE1, [y(x)], x = -5 .. 5, y = -5 .. 5, color = blue, scaling = constrained, axes = boxed)

 

 

DEplot(ODE1, y(x), x = -5 .. 5, y = -5 .. 5, [[y(0) = 2], [y(0) = 4]], linecolor = red, color = blue, scaling = constrained, axes = boxed)

 

Problem 2 (EULER METHOD)

NULL

NULL

ODE2 := diff(y(x), x) = sin(x*y(x))

diff(y(x), x) = sin(x*y(x))

(5)

 

 

dsolve(ODE2, y(x))

NULL

Maple returned no output! That means Maple is unable to solve the equation.

NULL

If you are curious what steps Maple went through to  find a solution before failing
to do so, you can ask to see the steps using the command "infolevel". The levels
of information that you can request range from 0 to 5.

 

````

dsolve(ODE2, y(x))

soln := dsolve({ODE2, y(0) = 3}, y(x), numeric)

proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if 1 < nargs then error "invalid input: too many arguments" end if; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then _xout := evalf[_EnvDSNumericSaveDigits](x_rkf45) else _xout := evalf(x_rkf45) end if; _dat := Array(1..4, {(1) = proc (_xin) local _xout, _dtbl, _dat, _vmap, _x0, _y0, _val, _dig, _n, _ne, _nd, _nv, _pars, _ini, _par, _i, _j, _k, _src; option `Copyright (c) 2002 by Waterloo Maple Inc. All rights reserved.`; table( [( "complex" ) = false ] ) _xout := _xin; _pars := []; _dtbl := array( 1 .. 4, [( 1 ) = (array( 1 .. 28, [( 1 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 2 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 3 ) = ([0, 0, 0, Array(1..0, 1..2, {}, datatype = float[8], order = C_order)]), ( 4 ) = (Array(1..65, {(1) = 1, (2) = 1, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 1, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 0, (14) = 0, (15) = 0, (16) = 0, (17) = 0, (18) = 21, (19) = 30000, (20) = 0, (21) = 0, (22) = 1, (23) = 4, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 1, (34) = 0, (35) = 0, (36) = 0, (37) = 0, (38) = 0, (39) = 0, (40) = 0, (41) = 0, (42) = 0, (43) = 1, (44) = 0, (45) = 0, (46) = 0, (47) = 0, (48) = 0, (49) = 0, (50) = 50, (51) = 1, (52) = 0, (53) = 0, (54) = 0, (55) = 0, (56) = 0, (57) = 0, (58) = 0, (59) = 10000, (60) = 0, (61) = 1000, (62) = 0, (63) = 0, (64) = -1, (65) = 0}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = .0, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = .0, (6) = .16290418802201975, (7) = .0, (8) = 0.10e-5, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = 1.0, (14) = .0, (15) = .49999999999999, (16) = .0, (17) = 1.0, (18) = 1.0, (19) = .0, (20) = .0, (21) = 1.0, (22) = 1.0, (23) = .0, (24) = .0, (25) = 0.10e-14, (26) = .0, (27) = .0, (28) = .0}, datatype = float[8], order = C_order)), ( 6 ) = (Array(1..1, {(1) = 3.0}, datatype = float[8], order = C_order)), ( 7 ) = ([Array(1..4, 1..7, {(1, 1) = .0, (1, 2) = .203125, (1, 3) = .3046875, (1, 4) = .75, (1, 5) = .8125, (1, 6) = .40625, (1, 7) = .8125, (2, 1) = 0.6378173828125e-1, (2, 2) = .0, (2, 3) = .279296875, (2, 4) = .27237892150878906, (2, 5) = -0.9686851501464844e-1, (2, 6) = 0.1956939697265625e-1, (2, 7) = .5381584167480469, (3, 1) = 0.31890869140625e-1, (3, 2) = .0, (3, 3) = -.34375, (3, 4) = -.335235595703125, (3, 5) = .2296142578125, (3, 6) = .41748046875, (3, 7) = 11.480712890625, (4, 1) = 0.9710520505905151e-1, (4, 2) = .0, (4, 3) = .40350341796875, (4, 4) = 0.20297467708587646e-1, (4, 5) = -0.6054282188415527e-2, (4, 6) = -0.4770040512084961e-1, (4, 7) = .77858567237854}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = 1.0, (2, 1) = .25, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = 1.0, (3, 1) = .1875, (3, 2) = .5625, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = 2.0, (4, 1) = .23583984375, (4, 2) = -.87890625, (4, 3) = .890625, (4, 4) = .0, (4, 5) = .0, (4, 6) = .2681884765625, (5, 1) = .1272735595703125, (5, 2) = -.5009765625, (5, 3) = .44921875, (5, 4) = -0.128936767578125e-1, (5, 5) = .0, (5, 6) = 0.626220703125e-1, (6, 1) = -0.927734375e-1, (6, 2) = .626220703125, (6, 3) = -.4326171875, (6, 4) = .1418304443359375, (6, 5) = -0.861053466796875e-1, (6, 6) = .3131103515625}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .386, (3) = .21, (4) = .63, (5) = 1.0, (6) = 1.0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .25, (2) = -.1043, (3) = .1035, (4) = -0.362e-1, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 1.544, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = .9466785280815533, (3, 2) = .25570116989825814, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = 3.3148251870684886, (4, 2) = 2.896124015972123, (4, 3) = .9986419139977808, (4, 4) = .0, (4, 5) = .0, (5, 1) = 1.2212245092262748, (5, 2) = 6.019134481287752, (5, 3) = 12.537083329320874, (5, 4) = -.687886036105895, (5, 5) = .0, (6, 1) = 1.2212245092262748, (6, 2) = 6.019134481287752, (6, 3) = 12.537083329320874, (6, 4) = -.687886036105895, (6, 5) = 1.0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = -5.6688, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = -2.4300933568337584, (3, 2) = -.20635991570891224, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = -.10735290581452621, (4, 2) = -9.594562251021896, (4, 3) = -20.470286148096154, (4, 4) = .0, (4, 5) = .0, (5, 1) = 7.496443313968615, (5, 2) = -10.246804314641219, (5, 3) = -33.99990352819906, (5, 4) = 11.708908932061595, (5, 5) = .0, (6, 1) = 8.083246795922411, (6, 2) = -7.981132988062785, (6, 3) = -31.52159432874373, (6, 4) = 16.319305431231363, (6, 5) = -6.0588182388340535}, datatype = float[8], order = C_order), Array(1..3, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 10.126235083446911, (2, 2) = -7.487995877607633, (2, 3) = -34.800918615557414, (2, 4) = -7.9927717075687275, (2, 5) = 1.0251377232956207, (3, 1) = -.6762803392806898, (3, 2) = 6.087714651678606, (3, 3) = 16.43084320892463, (3, 4) = 24.767225114183653, (3, 5) = -6.5943891257167815}, datatype = float[8], order = C_order)]), ( 9 ) = ([Array(1..1, {(1) = .1}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = 0}, datatype = integer[8]), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = 0}, datatype = integer[8])]), ( 8 ) = ([Array(1..1, {(1) = 3.0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..1, {(1, 1) = .0, (2, 0) = .0, (2, 1) = .0, (3, 0) = .0, (3, 1) = .0, (4, 0) = .0, (4, 1) = .0, (5, 0) = .0, (5, 1) = .0, (6, 0) = .0, (6, 1) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = y(x)]`; YP[1] := sin(X*Y[1]); 0 end proc, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]), ( 13 ) = (), ( 12 ) = (), ( 15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 18 ) = ([]), ( 19 ) = (0), ( 16 ) = ([0, 0, 0, 0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) option `[Y[1] = y(x)]`; YP[1] := sin(X*Y[1]); 0 end proc, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 27 ) = (""), ( 26 ) = (Array(1..0, {})), ( 25 ) = (Array(1..0, {})), ( 24 ) = (0), ( 28 ) = (0)  ] ))  ] ); _y0 := Array(0..1, {(1) = 0.}); _vmap := array( 1 .. 1, [( 1 ) = (1)  ] ); _x0 := _dtbl[1][5][5]; _n := _dtbl[1][4][1]; _ne := _dtbl[1][4][3]; _nd := _dtbl[1][4][4]; _nv := _dtbl[1][4][16]; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then if _Env_smart_dsolve_numeric = true or _dtbl[1][4][10] = 1 then if _xout = "left" then if type(_dtbl[2], 'table') then return _dtbl[2][5][1] end if elif _xout = "right" then if type(_dtbl[3], 'table') then return _dtbl[3][5][1] end if end if end if; return _dtbl[1][5][5] elif _xout = "method" then return _dtbl[1][15] elif _xout = "storage" then return evalb(_dtbl[1][4][10] = 1) elif _xout = "leftdata" then if not type(_dtbl[2], 'array') then return NULL else return eval(_dtbl[2]) end if elif _xout = "rightdata" then if not type(_dtbl[3], 'array') then return NULL else return eval(_dtbl[3]) end if elif _xout = "enginedata" then return eval(_dtbl[1]) elif _xout = "enginereset" then _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); return NULL elif _xout = "initial" then return procname(_y0[0]) elif _xout = "laxtol" then return _dtbl[`if`(member(_dtbl[4], {2, 3}), _dtbl[4], 1)][5][18] elif _xout = "numfun" then return `if`(member(_dtbl[4], {2, 3}), _dtbl[_dtbl[4]][4][18], 0) elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return procname(_y0[0]), [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _dtbl[4] <> 2 and _dtbl[4] <> 3 or _x0-_dtbl[_dtbl[4]][5][1] = 0. then error "no information is available on last computed point" else _xout := _dtbl[_dtbl[4]][5][1] end if elif _xout = "function" then if _dtbl[1][4][33]-2. = 0 then return eval(_dtbl[1][10], 1) else return eval(_dtbl[1][10][1], 1) end if elif _xout = "map" then return copy(_vmap) elif type(_xin, `=`) and type(rhs(_xin), 'list') and member(lhs(_xin), {"initial", "parameters", "initial_and_parameters"}) then _ini, _par := [], []; if lhs(_xin) = "initial" then _ini := rhs(_xin) elif lhs(_xin) = "parameters" then _par := rhs(_xin) elif select(type, rhs(_xin), `=`) <> [] then _par, _ini := selectremove(type, rhs(_xin), `=`) elif nops(rhs(_xin)) < nops(_pars)+1 then error "insufficient data for specification of initial and parameters" else _par := rhs(_xin)[-nops(_pars) .. -1]; _ini := rhs(_xin)[1 .. -nops(_pars)-1] end if; _xout := lhs(_xout); _i := false; if _par <> [] then _i := `dsolve/numeric/process_parameters`(_n, _pars, _par, _y0) end if; if _ini <> [] then _i := `dsolve/numeric/process_initial`(_n-_ne, _ini, _y0, _pars, _vmap) or _i end if; if _i then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars); if _Env_smart_dsolve_numeric = true and type(_y0[0], 'numeric') and _dtbl[1][4][10] <> 1 then procname("right") := _y0[0]; procname("left") := _y0[0] end if end if; if _xout = "initial" then return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)] elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] else return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)], [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] end if elif _xin = "eventstop" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then return 0 end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 <= _dtbl[5-_i][4][9] then _i := 5-_i; _dtbl[4] := _i; _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) elif 100 <= _dtbl[_i][4][9] then _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) else return 0 end if elif _xin = "eventstatus" then if _nv = 0 then error "this solution has no events" end if; _i := [selectremove(proc (a) options operator, arrow; _dtbl[1][3][1][a, 7] = 1 end proc, {seq(_j, _j = 1 .. round(_dtbl[1][3][1][_nv+1, 1]))})]; return ':-enabled' = _i[1], ':-disabled' = _i[2] elif _xin = "eventclear" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then error "no events to clear" end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 < _dtbl[5-_i][4][9] then _dtbl[4] := 5-_i; _i := 5-_i end if; if _dtbl[_i][4][9] < 100 then error "no events to clear" elif _nv < _dtbl[_i][4][9]-100 then error "event error condition cannot be cleared" else _j := _dtbl[_i][4][9]-100; if irem(round(_dtbl[_i][3][1][_j, 4]), 2) = 1 then error "retriggerable events cannot be cleared" end if; _j := round(_dtbl[_i][3][1][_j, 1]); for _k to _nv do if _dtbl[_i][3][1][_k, 1] = _j then if _dtbl[_i][3][1][_k, 2] = 3 then error "range events cannot be cleared" end if; _dtbl[_i][3][1][_k, 8] := _dtbl[_i][3][1][_nv+1, 8] end if end do; _dtbl[_i][4][17] := 0; _dtbl[_i][4][9] := 0; if _dtbl[1][4][10] = 1 then if _i = 2 then try procname(procname("left")) catch:  end try else try procname(procname("right")) catch:  end try end if end if end if; return  elif type(_xin, `=`) and member(lhs(_xin), {"eventdisable", "eventenable"}) then if _nv = 0 then error "this solution has no events" end if; if type(rhs(_xin), {('list')('posint'), ('set')('posint')}) then _i := {op(rhs(_xin))} elif type(rhs(_xin), 'posint') then _i := {rhs(_xin)} else error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; if select(proc (a) options operator, arrow; _nv < a end proc, _i) <> {} then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _k := {}; for _j to _nv do if member(round(_dtbl[1][3][1][_j, 1]), _i) then _k := `union`(_k, {_j}) end if end do; _i := _k; if lhs(_xin) = "eventdisable" then _dtbl[4] := 0; _j := [evalb(assigned(_dtbl[2]) and member(_dtbl[2][4][17], _i)), evalb(assigned(_dtbl[3]) and member(_dtbl[3][4][17], _i))]; for _k in _i do _dtbl[1][3][1][_k, 7] := 0; if assigned(_dtbl[2]) then _dtbl[2][3][1][_k, 7] := 0 end if; if assigned(_dtbl[3]) then _dtbl[3][3][1][_k, 7] := 0 end if end do; if _j[1] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[2][3][4][_k, 1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to defined init `, _dtbl[2][3][4][_k, 1]); _dtbl[2][3][1][_k, 8] := _dtbl[2][3][4][_k, 1] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to rate hysteresis init `, _dtbl[2][5][24]); _dtbl[2][3][1][_k, 8] := _dtbl[2][5][24] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to initial init `, _x0); _dtbl[2][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to fireinitial init `, _x0-1); _dtbl[2][3][1][_k, 8] := _x0-1 end if end do; _dtbl[2][4][17] := 0; _dtbl[2][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("left")) end if end if; if _j[2] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[3][3][4][_k, 2], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to defined init `, _dtbl[3][3][4][_k, 2]); _dtbl[3][3][1][_k, 8] := _dtbl[3][3][4][_k, 2] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to rate hysteresis init `, _dtbl[3][5][24]); _dtbl[3][3][1][_k, 8] := _dtbl[3][5][24] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to initial init `, _x0); _dtbl[3][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to fireinitial init `, _x0+1); _dtbl[3][3][1][_k, 8] := _x0+1 end if end do; _dtbl[3][4][17] := 0; _dtbl[3][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("right")) end if end if else for _k in _i do _dtbl[1][3][1][_k, 7] := 1 end do; _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); _dtbl[4] := 0; if _dtbl[1][4][10] = 1 then if _x0 <= procname("right") then try procname(procname("right")) catch:  end try end if; if procname("left") <= _x0 then try procname(procname("left")) catch:  end try end if end if end if; return  elif type(_xin, `=`) and lhs(_xin) = "eventfired" then if not type(rhs(_xin), 'list') then error "'eventfired' must be specified as a list" end if; if _nv = 0 then error "this solution has no events" end if; if _dtbl[4] <> 2 and _dtbl[4] <> 3 then error "'direction' must be set prior to calling/setting 'eventfired'" end if; _i := _dtbl[4]; _val := NULL; if not assigned(_EnvEventRetriggerWarned) then _EnvEventRetriggerWarned := false end if; for _k in rhs(_xin) do if type(_k, 'integer') then _src := _k elif type(_k, 'integer' = 'anything') and type(evalf(rhs(_k)), 'numeric') then _k := lhs(_k) = evalf[max(Digits, 18)](rhs(_k)); _src := lhs(_k) else error "'eventfired' entry is not valid: %1", _k end if; if _src < 1 or round(_dtbl[1][3][1][_nv+1, 1]) < _src then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _src := {seq(`if`(_dtbl[1][3][1][_j, 1]-_src = 0., _j, NULL), _j = 1 .. _nv)}; if nops(_src) <> 1 then error "'eventfired' can only be set/queried for root-finding events and time/interval events" end if; _src := _src[1]; if _dtbl[1][3][1][_src, 2] <> 0. and _dtbl[1][3][1][_src, 2]-2. <> 0. then error "'eventfired' can only be set/queried for root-finding events and time/interval events" elif irem(round(_dtbl[1][3][1][_src, 4]), 2) = 1 then if _EnvEventRetriggerWarned = false then WARNING(`'eventfired' has no effect on events that retrigger`) end if; _EnvEventRetriggerWarned := true end if; if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then _val := _val, undefined elif type(_dtbl[_i][3][4][_src, _i-1], 'undefined') or _i = 2 and _dtbl[2][3][1][_src, 8] < _dtbl[2][3][4][_src, 1] or _i = 3 and _dtbl[3][3][4][_src, 2] < _dtbl[3][3][1][_src, 8] then _val := _val, _dtbl[_i][3][1][_src, 8] else _val := _val, _dtbl[_i][3][4][_src, _i-1] end if; if type(_k, `=`) then if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then error "cannot set event code for a rate hysteresis event" end if; userinfo(3, {'events', 'eventreset'}, `manual set event code `, _src, ` to value `, rhs(_k)); _dtbl[_i][3][1][_src, 8] := rhs(_k); _dtbl[_i][3][4][_src, _i-1] := rhs(_k) end if end do; return [_val] elif type(_xin, `=`) and lhs(_xin) = "direction" then if not member(rhs(_xin), {-1, 1, ':-left', ':-right'}) then error "'direction' must be specified as either '1' or 'right' (positive) or '-1' or 'left' (negative)" end if; _src := `if`(_dtbl[4] = 2, -1, `if`(_dtbl[4] = 3, 1, undefined)); _i := `if`(member(rhs(_xin), {1, ':-right'}), 3, 2); _dtbl[4] := _i; _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if; return _src elif _xin = "eventcount" then if _dtbl[1][3][1] = 0 or _dtbl[4] <> 2 and _dtbl[4] <> 3 then return 0 else return round(_dtbl[_dtbl[4]][3][1][_nv+1, 12]) end if elif type(_xin, `=`) and lhs(_xin) = "setdatacallback" then if not type(rhs(_xin), 'nonegint') then error "data callback must be a nonnegative integer (address)" end if; _dtbl[1][28] := rhs(_xin) else return "procname" end if end if; if _xout = _x0 then return [_x0, seq(evalf(_dtbl[1][6][_vmap[_i]]), _i = 1 .. _n-_ne)] end if; _i := `if`(_x0 <= _xout, 3, 2); if _xin = "last" and 0 < _dtbl[_i][4][9] and _dtbl[_i][4][9] < 100 then _dat := eval(_dtbl[_i], 2); _j := _dat[4][20]; return [_dat[11][_j, 0], seq(_dat[11][_j, _vmap[_i]], _i = 1 .. _n-_ne-_nd), seq(_dat[8][1][_vmap[_i]], _i = _n-_ne-_nd+1 .. _n-_ne)] end if; if not type(_dtbl[_i], 'array') then _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if end if; if _xin <> "last" then if 0 < 0 then if `dsolve/numeric/checkglobals`(op(_dtbl[1][14]), _pars, _n, _y0) then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars, _i) end if end if; if _dtbl[1][4][7] = 0 then error "parameters must be initialized before solution can be computed" end if end if; _dat := eval(_dtbl[_i], 2); _dtbl[4] := _i; try _src := `dsolve/numeric/SC/IVPrun`(_dat, _xout) catch: userinfo(2, `dsolve/debug`, print(`Exception in solnproc:`, [lastexception][2 .. -1])); error  end try; if _dat[17] <> _dtbl[1][17] then _dtbl[1][17] := _dat[17]; _dtbl[1][10] := _dat[10] end if; if _src = 0 and 100 < _dat[4][9] then _val := _dat[3][1][_nv+1, 8] else _val := _dat[11][_dat[4][20], 0] end if; if _src <> 0 or _dat[4][9] <= 0 then _dtbl[1][5][1] := _xout else _dtbl[1][5][1] := _val end if; if _i = 3 and _val < _xout then Rounding := -infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further right of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further right of %1, maxfun limit exceeded (see <a href='http://www.maplesoft.com/support/help/search.aspx?term=dsolve,maxfun' target='_new'>?dsolve,maxfun</a> for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further right of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further right of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further right of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further right of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further right of %1", evalf[8](_val) end if elif _i = 2 and _xout < _val then Rounding := infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further left of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further left of %1, maxfun limit exceeded (see <a href='http://www.maplesoft.com/support/help/search.aspx?term=dsolve,maxfun' target='_new'>?dsolve,maxfun</a> for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further left of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further left of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further left of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further left of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further left of %1", evalf[8](_val) end if end if; if _EnvInFsolve = true then _dig := _dat[4][26]; if type(_EnvDSNumericSaveDigits, 'posint') then _dat[4][26] := _EnvDSNumericSaveDigits else _dat[4][26] := Digits end if; _Env_dsolve_SC_native := true; if _dat[4][25] = 1 then _i := 1; _dat[4][25] := 2 else _i := _dat[4][25] end if; _val := `dsolve/numeric/SC/IVPval`(_dat, _xout, _src); _dat[4][25] := _i; _dat[4][26] := _dig; [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] else Digits := _dat[4][26]; _val := `dsolve/numeric/SC/IVPval`(eval(_dat, 2), _xout, _src); [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] end if end proc, (2) = Array(0..0, {}), (3) = [x, y(x)], (4) = []}); _vars := _dat[3]; _pars := map(lhs, _dat[4]); _n := nops(_vars)-1; _solnproc := _dat[1]; if not type(_xout, 'numeric') then if member(x_rkf45, ["start", 'start', "method", 'method', "left", 'left', "right", 'right', "leftdata", "rightdata", "enginedata", "eventstop", 'eventstop', "eventclear", 'eventclear', "eventstatus", 'eventstatus', "eventcount", 'eventcount', "laxtol", 'laxtol', "numfun", 'numfun', NULL]) then _res := _solnproc(convert(x_rkf45, 'string')); if 1 < nops([_res]) then return _res elif type(_res, 'array') then return eval(_res, 1) elif _res <> "procname" then return _res end if elif member(x_rkf45, ["last", 'last', "initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x_rkf45, 'string'); _res := _solnproc(_xout); if _xout = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] end if elif type(_xout, `=`) and member(lhs(_xout), ["initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(lhs(x_rkf45), 'string') = rhs(x_rkf45); if type(rhs(_xout), 'list') then _res := _solnproc(_xout) else error "initial and/or parameter values must be specified in a list" end if; if lhs(_xout) = "initial" then return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] elif lhs(_xout) = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["eventdisable", 'eventdisable', "eventenable", 'eventenable', "eventfired", 'eventfired', "direction", 'direction', NULL]) then return _solnproc(convert(lhs(x_rkf45), 'string') = rhs(x_rkf45)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _vars end if; if procname <> unknown then return ('procname')(x_rkf45) else _ndsol := 1; _ndsol := _ndsol; _ndsol := pointto(_dat[2][0]); return ('_ndsol')(x_rkf45) end if end if; try _res := _solnproc(_xout); [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] catch: error  end try end proc

(6)

[soln(1), soln(1.5), soln(2), soln(2.5), soln(3), soln(3.5), soln(4), soln(4.5), soln(5)]

[[x = 1., y(x) = HFloat(3.5280317971877286)], [x = 1.5, y(x) = HFloat(3.1226400064202693)], [x = 2., y(x) = HFloat(2.653970420106217)], [x = 2.5, y(x) = HFloat(2.323175669304077)], [x = 3., y(x) = HFloat(2.3019187052877816)], [x = 3.5, y(x) = HFloat(2.6587033583656714)], [x = 4., y(x) = HFloat(2.497876146260152)], [x = 4.5, y(x) = HFloat(2.220305976552359)], [x = 5., y(x) = HFloat(1.9758297601444128)]]

(7)

p1 := odeplot(soln, [x, y(x)], x = 0 .. 5, color = magenta, thickness = 2, scaling = constrained, view = [0 .. 5, 0 .. 5])

 

p2 := dfieldplot(ODE2, [y(x)], x = 0 .. 5, y = 0 .. 5, color = blue, scaling = constrained)

display(p1, p2)

 

DEplot(ODE2, y(x), x = 0 .. 5, y = 0 .. 5, [[y(0) = 3]], linecolor = red, color = blue, scaling = constrained, axes = boxed)

 

Construct approximate solutions for x from 0 to 5 to the initial problem 2 using Euler's method with the three different step sizes Δx=0.5, 0.25, 0.125

f := proc (x, y) options operator, arrow; sin(x*y) end proc

proc (x, y) options operator, arrow; sin(y*x) end proc

(8)

Eulermethod := proc (f, x_start, y_start, dx, n_total) local x, y, Y, i, k; x[1] := x_start; y[1] := y_start; for i to n_total do y[i+1] := y[i]+f(x[i], y[i])*dx; x[i+1] := x[1]+i*dx end do; Y := [seq([x[k], y[k]], k = 1 .. n_total)]; return Y end proc

proc (f, x_start, y_start, dx, n_total) local x, y, Y, i, k; x[1] := x_start; y[1] := y_start; for i to n_total do y[i+1] := y[i]+f(x[i], y[i])*dx; x[i+1] := x[1]+i*dx end do; Y := [seq([x[k], y[k]], k = 1 .. n_total)]; return Y end proc

(9)

NULL

NULL

NULL

NULL

A1 := Eulermethod(f, 0, 3, .5, 10)

[[0, 3], [.5, 3.], [1.0, 3.498747493], [1.5, 3.323942476], [2.0, 2.842530099], [2.5, 2.560983065], [3.0, 2.620477947], [3.5, 3.120464063], [4.0, 2.621830593], [4.5, 2.185032319]]

(10)

A2 := Eulermethod(f, 0, 3, .25, 20)

[[0, 3], [.25, 3.], [.50, 3.170409690], [.75, 3.420383740], [1.00, 3.556616077], [1.25, 3.455813236], [1.50, 3.224836021], [1.75, 2.976782400], [2.00, 2.757025820], [2.25, 2.583147563], [2.50, 2.469680146], [2.75, 2.442487816], [3.00, 2.547535655], [3.25, 2.791971512], [3.50, 2.877900373], [3.75, 2.727027361], [4.00, 2.547414295], [4.25, 2.374301829], [4.50, 2.219839448], [4.75, 2.086091178]]

(11)

A3 := Eulermethod(f, 0, 3, .125, 40)

[[0, 3], [.125, 3.], [.250, 3.045784066], [.375, 3.132030172], [.500, 3.247342837], [.625, 3.372168142], [.750, 3.479586272], [.875, 3.542983060], [1.000, 3.548166882], [1.125, 3.498733738], [1.250, 3.409546073], [1.375, 3.297015011], [1.500, 3.174012084], [1.625, 3.049159854], [1.750, 2.927817142], [1.875, 2.813241461], [2.000, 2.707496816], [2.125, 2.612101612], [2.250, 2.528513147], [2.375, 2.458549923], [2.500, 2.404840953], [2.625, 2.371369082], [2.750, 2.364080535], [2.875, 2.391119623], [3.000, 2.460798019], [3.125, 2.572154041], [3.250, 2.695044010], [3.375, 2.772263417], [3.500, 2.780805371], [3.625, 2.742906337], [3.750, 2.680985435], [3.875, 2.607451728], [4.000, 2.528940353], [4.125, 2.449278434], [4.250, 2.370825619], [4.375, 2.295054886], [4.500, 2.222824117], [4.625, 2.154537647], [4.750, 2.090275081], [4.875, 2.029905439]]

(12)

NULL

p3 := pointplot(A1, color = green, scaling = constrained, symbol = circle)

p4 := pointplot(A2, color = black, scaling = constrained, symbol = asterisk)

p5 := pointplot(A3, color = red, scaling = constrained, symbol = diamond)

NULL

display([p1, p3, p4, p5])

 

NULL

In this plot, the differences between the lines visually represent how using different step sizes affects the overall solution accuracy. The red points are closest to what a more accurate numerical solution would look like, while the green points are more of a rough approximation.

 

 

Problem 3 (IMPROVED EULER METHOD

``

 

 

 

ImprovedEulermethod := proc (f, x_start, y_start, dx, n_total) local x, y, Y, i, k, k1, k2; x[1] := x_start; y[1] := y_start; for i to n_total do k1 := f(x[i], y[i]); k2 := f(x[i]+dx, y[i]+k1*dx); y[i+1] := y[i]+(1/2)*dx*(k1+k2); x[i+1] := x[1]+i*dx end do; Y := [seq([x[k], y[k]], k = 1 .. n_total+1)]; return Y end proc

proc (f, x_start, y_start, dx, n_total) local x, y, Y, i, k, k1, k2; x[1] := x_start; y[1] := y_start; for i to n_total do k1 := f(x[i], y[i]); k2 := f(x[i]+dx, y[i]+k1*dx); y[i+1] := y[i]+(1/2)*dx*(k1+k2); x[i+1] := x[1]+i*dx end do; Y := [seq([x[k], y[k]], k = 1 .. n_total+1)]; return Y end proc

(13)

A := ImprovedEulermethod(f, 0, 3, .5, 10); B := ImprovedEulermethod(f, 0, 3, .25, 20); C := ImprovedEulermethod(f, 0, 3, .125, 40)

[[0, 3], [.125, 3.022892033], [.250, 3.089335383], [.375, 3.190999573], [.500, 3.311460714], [.625, 3.426126738], [.750, 3.508312296], [.875, 3.539992283], [1.000, 3.518184043], [1.125, 3.451931748], [1.250, 3.354946855], [1.375, 3.240089444], [1.500, 3.117181611], [1.625, 2.992925708], [1.750, 2.871597591], [1.875, 2.755823574], [2.000, 2.647215450], [2.125, 2.546845751], [2.250, 2.455612702], [2.375, 2.374560360], [2.500, 2.305227222], [2.625, 2.250113264], [2.750, 2.213376654], [2.875, 2.201814459], [3.000, 2.225589190], [3.125, 2.295322044], [3.250, 2.406184220], [3.375, 2.516909932], [3.500, 2.583378006], [3.625, 2.599915321], [3.750, 2.579967129], [3.875, 2.537160528], [4.000, 2.481052245], [4.125, 2.417781798], [4.250, 2.351265142], [4.375, 2.284028571], [4.500, 2.217702881], [4.625, 2.153313287], [4.750, 2.091461577], [4.875, 2.032452273], [5.000, 1.976386380]]

(14)

NULL

plot2 := pointplot(A, color = green, scaling = constrained, symbol = circle)

plot3 := pointplot(B, color = black, scaling = constrained, symbol = asterisk)

plot4 := pointplot(C, color = red, scaling = constrained, symbol = diamond)

NULL

display([p1, plot2, plot3, plot4])

 

The Improved Euler method, by averaging slopes, provides a significant improvement over the basic Euler method, particularly when the step size is relatively large.


Download Diff_equations.mw

In this post, I would like to share some exercises that I recently taught to an undergraduate student using Maple. These exercises aimed to deepen their understanding of mathematical concepts through computational exploration and visualization. With its powerful symbolic computation capabilities, Maple proved to be an excellent tool for this purpose. Below, I present a few of the exercises and the insights they provided. Interestingly, the student found Maple to be more user-friendly and efficient compared to the software he usually uses for his studies. Below, I present a few of the exercises and the insights they provided.

One of the first topics we tackled was the Fourier series. We used Maple to illustrate how the Fourier series approximates a given function as more terms are added. We explored this through both static plots and interactive animations.

To help the student understand the behavior of different types of functions, we defined piecewise functions using Maple's piecewise command. This allowed us to model functions that behave differently over various intervals, such as the following cubic function exercise

Maple's Explore command was an effective tool for creating an interactive learning environment. We used it to create sliders that allowed the student to vary parameters, such as the number of terms in a Fourier series, and see the immediate impact on the plot.

restart; with(plots)

" F(x):={[[-1,-1<x<0],[1,0<x<1]];  "

proc (x) options operator, arrow, function_assign; piecewise(-1 < x and x < 0, -1, 0 < x and x < 1, 1) end proc

(1)

p1 := plot(piecewise(-3 < x and x < -1, F(x+2), -1 < x and x < 1, F(x), 1 < x and x < 3, F(x-2)), x = -3 .. 3, color = blue)

 

L := 2; a__0 := (int(F(x), x = -(1/2)*L .. (1/2)*L))/L

0

(2)

a__n := proc (n) options operator, arrow; 2*(int(F(x)*cos(2*n*Pi*x/L), x = -(1/2)*L .. (1/2)*L))/L end proc

proc (n) options operator, arrow; 2*(int(F(x)*cos(2*n*Pi*x/L), x = -(1/2)*L .. (1/2)*L))/L end proc

(3)

b__n := proc (n) options operator, arrow; 2*(int(F(x)*sin(2*n*Pi*x/L), x = -(1/2)*L .. (1/2)*L))/L end proc

proc (n) options operator, arrow; 2*(int(F(x)*sin(2*n*Pi*x/L), x = -(1/2)*L .. (1/2)*L))/L end proc

(4)

" #` Fourier series function`  `F__fourier`(x,N):=`a__0`+(&sum;)(`a__n`(n)&lowast;cos(2 * n * Pi * x / L) +`b__n`(n)&lowast;sin(2* n * Pi * x / L));"

proc (x, N) options operator, arrow, function_assign; a__0+sum(a__n(n)*cos(2*n*Pi*x/L)+b__n(n)*sin(2*n*Pi*x/L), n = 1 .. N) end proc

(5)

p2 := plot([F__fourier(x, 40)], x = -3 .. 3, numpoints = 200, color = [purple])

display([p1, p2])

 

Explore(plot([piecewise(-3 < x and x < -1, F(x+2), -1 < x and x < 1, F(x), 1 < x and x < 3, F(x-2)), F__fourier(x, N)], x = -3 .. 3, color = [blue, purple], numpoints = 200), N = 1 .. 40, title = "Fourier Series Approximation with N Terms")

restart; with(plots)

" #` Define the piecewise function`  F(x):={[[0,-1<x<0],[x^(2),0<x<1]];  "

proc (x) options operator, arrow, function_assign; piecewise(-1 < x and x < 0, 0, 0 < x and x < 1, x^2) end proc

(6)

p3 := plot(piecewise(-3 < x and x < -1, F(x+2), -1 < x and x < 1, F(x), 1 < x and x < 3, F(x-2)), x = -3 .. 3, color = blue)

 

L := 2; a__0 := (int(F(x), x = -(1/2)*L .. (1/2)*L))/L

1/6

(7)

a__n := proc (n) options operator, arrow; 2*(int(F(x)*cos(2*n*Pi*x/L), x = -(1/2)*L .. (1/2)*L))/L end proc

proc (n) options operator, arrow; 2*(int(F(x)*cos(2*n*Pi*x/L), x = -(1/2)*L .. (1/2)*L))/L end proc

(8)

b__n := proc (n) options operator, arrow; 2*(int(F(x)*sin(2*n*Pi*x/L), x = -(1/2)*L .. (1/2)*L))/L end proc

proc (n) options operator, arrow; 2*(int(F(x)*sin(2*n*Pi*x/L), x = -(1/2)*L .. (1/2)*L))/L end proc

(9)

" #` Fourier series function`  `F__fourier`(x,N):=`a__0`+(&sum;)(`a__n`(n)&lowast;cos(2 * n * Pi * x / L) +`b__n`(n)&lowast;sin(2* n * Pi * x / L));"

proc (x, N) options operator, arrow, function_assign; a__0+sum(a__n(n)*cos(2*n*Pi*x/L)+b__n(n)*sin(2*n*Pi*x/L), n = 1 .. N) end proc

(10)

p4 := plot([F__fourier(x, 40)], x = -3 .. 3, numpoints = 200, color = [purple])

display([p3, p4])

 

Explore(plot([piecewise(-3 < x and x < -1, F(x+2), -1 < x and x < 1, F(x), 1 < x and x < 3, F(x-2)), F__fourier(x, N)], x = -3 .. 3, color = [blue, purple], numpoints = 200), N = 1 .. 40, title = "Fourier Series Approximation with N Terms")

restart; with(plots)

" #` Define the piecewise function`  F(x):={[[x+2,-2<x<0],[2-2 x,0<x<2]]; "

proc (x) options operator, arrow, function_assign; piecewise(-2 < x and x < 0, x+2, 0 < x and x < 2, 2-2*x) end proc

(11)

p5 := plot(piecewise(-3 < x and x < -1, F(x+2), -1 < x and x < 1, F(x), 1 < x and x < 3, F(x-2)), x = -3 .. 3, color = blue)

 

L := 2; a__0 := (int(F(x), x = -(1/2)*L .. (1/2)*L))/L

5/4

(12)

a__n := proc (n) options operator, arrow; 2*(int(F(x)*cos(2*n*Pi*x/L), x = -(1/2)*L .. (1/2)*L))/L end proc

proc (n) options operator, arrow; 2*(int(F(x)*cos(2*n*Pi*x/L), x = -(1/2)*L .. (1/2)*L))/L end proc

(13)

b__n := proc (n) options operator, arrow; 2*(int(F(x)*sin(2*n*Pi*x/L), x = -(1/2)*L .. (1/2)*L))/L end proc

proc (n) options operator, arrow; 2*(int(F(x)*sin(2*n*Pi*x/L), x = -(1/2)*L .. (1/2)*L))/L end proc

(14)

" #` Fourier series function`  `F__fourier`(x,N):=`a__0`+(&sum;)(`a__n`(n)&lowast;cos(2 * n * Pi * x / L) +`b__n`(n)&lowast;sin(2* n * Pi * x / L));"

proc (x, N) options operator, arrow, function_assign; a__0+sum(a__n(n)*cos(2*n*Pi*x/L)+b__n(n)*sin(2*n*Pi*x/L), n = 1 .. N) end proc

(15)

p6 := plot([F__fourier(x, 40)], x = -3 .. 3, numpoints = 200, color = [purple])

display([p5, p6])

 

Explore(plot([piecewise(-3 < x and x < -1, F(x+2), -1 < x and x < 1, F(x), 1 < x and x < 3, F(x-2)), F__fourier(x, N)], x = -3 .. 3, color = [blue, purple], numpoints = 200), N = 1 .. 40, title = "Fourier Series Approximation with N Terms")

restart; with(plots)

F := proc (x) options operator, arrow; piecewise(-1 < x and x < 1, x-x^3, 0) end proc

proc (x) options operator, arrow; piecewise(-1 < x and x < 1, x-x^3, 0) end proc

(16)

p7 := plot(piecewise(-3 < x and x < -1, F(x+2), -1 < x and x < 1, F(x), 1 < x and x < 3, F(x-2)), x = -3 .. 3, color = blue)

 

L := 2; a__0 := (int(F(x), x = -(1/2)*L .. (1/2)*L))/L

0

(17)

a__n := proc (n) options operator, arrow; 2*(int(F(x)*cos(2*n*Pi*x/L), x = -(1/2)*L .. (1/2)*L))/L end proc

proc (n) options operator, arrow; 2*(int(F(x)*cos(2*n*Pi*x/L), x = -(1/2)*L .. (1/2)*L))/L end proc

(18)

b__n := proc (n) options operator, arrow; 2*(int(F(x)*sin(2*n*Pi*x/L), x = -(1/2)*L .. (1/2)*L))/L end proc

proc (n) options operator, arrow; 2*(int(F(x)*sin(2*n*Pi*x/L), x = -(1/2)*L .. (1/2)*L))/L end proc

(19)

b__n(n)

-4*(n^2*Pi^2*sin(n*Pi)+3*cos(n*Pi)*Pi*n-3*sin(n*Pi))/(n^4*Pi^4)

(20)

" #` Fourier series function`  `F__fourier`(x,N):=`a__0`+(&sum;)(`a__n`(n)&lowast;cos(2 * n * Pi * x / L) +`b__n`(n)&lowast;sin(2* n * Pi * x / L));      #` Plot the Fourier series approximation`  p8:=plot([`F__fourier`(x,40)],x = -3.. 3 ,numpoints=200, color=[red]) :"

proc (x, N) options operator, arrow, function_assign; a__0+sum(a__n(n)*cos(2*n*Pi*x/L)+b__n(n)*sin(2*n*Pi*x/L), n = 1 .. N) end proc

(21)

display([p7, p8])

 

Explore(plot([piecewise(-3 < x and x < -1, F(x+2), -1 < x and x < 1, F(x), 1 < x and x < 3, F(x-2)), F__fourier(x, N)], x = -3 .. 3, color = [blue, purple], numpoints = 200), N = 1 .. 40, title = "Fourier Series Approximation with N Terms")

NULL

Download Fourier_Series.mw

Today in class, we presented an exercise based on the paper titled "Analysis of regular and chaotic dynamics in a stochastic eco-epidemiological model" by Bashkirtseva, Ryashko, and Ryazanova (2020). In this exercise, we kept all parameters of the model the same as in the paper, but we varied the parameter β, which represents the rate of infection spread. The goal was to observe how changes in β impact the system's dynamics, particularly focusing on the transition between regular and chaotic behavior.

This exercise involves studying a mathematical model that appears in eco-epidemiology. The model is described by the following set of equations:

dx/dt = rx-bx^2-cxy-`&beta;xz`/(a+x)-a[1]*yz/(e+y)

dy/dt = -`&mu;y`+`&beta;xy`/(a+x)-a[2]*yz/(d+y)

" (dz)/(dt)=-mz+((c[1 ]a[1])[ ]xz)/(e[]+x)+((c[1 ]a[2])[ ]yz)/(d+y)"

 

where r, b, c, β, α,a[1],a[2], e, d, m, c[1], c[2], μ>0 are given parameters. This model generalizes the classic predator-prey system by incorporating disease dynamics within the prey population. The populations are divided into the following groups:

 

• 

Susceptible prey population (x): Individuals in the prey population that are healthy but can become infected by a disease.

• 

Infected prey population (y): Individuals in the prey population that are infected and can transmit the disease to others.

• 

Predator population (z): The predator population that feeds on both susceptible (x) and infected (y) prey.

 

The initial conditions are always x(0)=0.2, y(0)=0.05, z(0)=0.05,  and we will vary the parameter β.;

 

For this exercise, the parameters are fixed as follows:

 

"r=1,` b`=1,` c`=0.01, a=0.36 ,` a`[1]=0.01,` a`[2]=0.05,` e`[]=15,` m`=0.01,` d`=0.5,` c`[1]=2,` `c[2]==1,` mu`=0.4."

NULL

Task (a)

• 

Solve the system numerically for the given parameter values and initial conditions with β=0.6 over the time interval t2[0,20000].

• 

Plot the solutions x(t), y(t), and  z(t) over this time interval.

• 

Comment on the model's predictions, keeping in mind that the time units are usually days.

• 

Also, plot the trajectory in the 3D space (x,y,z).

 

restart

r := 1; b := 1; f := 0.1e-1; alpha := .36; a[1] := 0.1e-1; a[2] := 0.5e-1; e := 15; m := 0.1e-1; d := .5; c[1] := 2; c[2] := 1; mu := .4; beta := .6

sys := {diff(x(t), t) = r*x(t)-b*x(t)^2-f*x(t)*y(t)-beta*x(t)*y(t)/(alpha+x(t))-a[1]*x(t)*z(t)/(e+x(t)), diff(y(t), t) = -mu*y(t)+beta*x(t)*y(t)/(alpha+x(t))-a[2]*y(t)*z(t)/(d+y(t)), diff(z(t), t) = -m*z(t)+c[1]*a[1]*x(t)*z(t)/(e+x(t))+c[2]*a[2]*y(t)*z(t)/(d+y(t))}

{diff(x(t), t) = x(t)-x(t)^2-0.1e-1*x(t)*y(t)-.6*x(t)*y(t)/(.36+x(t))-0.1e-1*x(t)*z(t)/(15+x(t)), diff(y(t), t) = -.4*y(t)+.6*x(t)*y(t)/(.36+x(t))-0.5e-1*y(t)*z(t)/(.5+y(t)), diff(z(t), t) = -0.1e-1*z(t)+0.2e-1*x(t)*z(t)/(15+x(t))+0.5e-1*y(t)*z(t)/(.5+y(t))}

(1)

ics := {x(0) = .2, y(0) = 0.5e-1, z(0) = 0.5e-1}

{x(0) = .2, y(0) = 0.5e-1, z(0) = 0.5e-1}

(2)

NULL

sol := dsolve(`union`(sys, ics), {x(t), y(t), z(t)}, numeric, range = 0 .. 20000, maxfun = 0, output = listprocedure, abserr = 0.1e-7, relerr = 0.1e-7)

`[Length of output exceeds limit of 1000000]`

(3)

X := subs(sol, x(t)); Y := subs(sol, y(t)); Z := subs(sol, z(t))

``

plot('[X(t)]', t = 0 .. 20000, numpoints = 350, title = "Trajectory of x(t)", axes = boxed, gridlines, color = ["#40e0d0"])

 

plot('[Y(t)]', t = 0 .. 20000, numpoints = 350, title = "Trajectory", axes = boxed, gridlines, title = "Trajectory of y(t)", color = ["SteelBlue"])

 

``

plot('[Z(t)]', t = 0 .. 20000, numpoints = 350, title = "Trajectory", axes = boxed, gridlines, title = "Trajectory of Z(t)", color = "Black"); with(DEtools)

 

with(DEtools)

DEplot3d(sys, {x(t), y(t), z(t)}, t = 0 .. 20000, [[x(0) = .2, y(0) = 0.5e-1, z(0) = 0.5e-1]], numpoints = 35000, color = blue, thickness = 1, linestyle = solid)

 

Task (b)

• 

Repeat the study in part (a) with the same initial conditions but set β=0.61.

NULL

restart

r := 1; b := 1; f := 0.1e-1; alpha := .36; a[1] := 0.1e-1; a[2] := 0.5e-1; e := 15; m := 0.1e-1; d := .5; c[1] := 2; c[2] := 1; mu := .4; beta := .61

NULL

sys := {diff(x(t), t) = r*x(t)-b*x(t)^2-f*x(t)*y(t)-beta*x(t)*y(t)/(alpha+x(t))-a[1]*x(t)*z(t)/(e+x(t)), diff(y(t), t) = -mu*y(t)+beta*x(t)*y(t)/(alpha+x(t))-a[2]*y(t)*z(t)/(d+y(t)), diff(z(t), t) = -m*z(t)+c[1]*a[1]*x(t)*z(t)/(e+x(t))+c[2]*a[2]*y(t)*z(t)/(d+y(t))}

{diff(x(t), t) = x(t)-x(t)^2-0.1e-1*x(t)*y(t)-.61*x(t)*y(t)/(.36+x(t))-0.1e-1*x(t)*z(t)/(15+x(t)), diff(y(t), t) = -.4*y(t)+.61*x(t)*y(t)/(.36+x(t))-0.5e-1*y(t)*z(t)/(.5+y(t)), diff(z(t), t) = -0.1e-1*z(t)+0.2e-1*x(t)*z(t)/(15+x(t))+0.5e-1*y(t)*z(t)/(.5+y(t))}

(4)

NULL

ics := {x(0) = .2, y(0) = 0.5e-1, z(0) = 0.5e-1}

{x(0) = .2, y(0) = 0.5e-1, z(0) = 0.5e-1}

(5)

sol := dsolve(`union`(sys, ics), {x(t), y(t), z(t)}, numeric, range = 0 .. 20000, maxfun = 0, output = listprocedure, abserr = 0.1e-7, relerr = 0.1e-7)

`[Length of output exceeds limit of 1000000]`

(6)

X := subs(sol, x(t)); Y := subs(sol, y(t)); Z := subs(sol, z(t))

NULL

plot('[X(t)]', t = 0 .. 20000, numpoints = 350, title = "Trajectory of x(t)", axes = boxed, gridlines, color = ["Blue"])

 

plot('[Y(t)]', t = 0 .. 20000, numpoints = 350, title = "Trajectory of  Y(t)", axes = boxed, gridlines, color = "Red")

 

plot('[Z(t)]', t = 0 .. 20000, numpoints = 350, title = "Trajectory of  Y(t)", axes = boxed, gridlines, color = "Black")

 

NULL

with(DEtools)

DEplot3d(sys, {x(t), y(t), z(t)}, t = 0 .. 20000, [[x(0) = .2, y(0) = 0.5e-1, z(0) = 0.5e-1]], maxfun = 0, numpoints = 35000, color = blue, thickness = 1, linestyle = solid)

 

The rate of the infection spread is affected by the average number of contacts each person has (β=0.6) and increases depending on the degree of transmission within the population, in particular within specific subpopulations (such as those in rural areas). A detailed epidemiological study showed that the spread of infection is most significant in urban areas, where population density is higher, while in rural areas, the rate of infection remains relatively low. This suggests that additional public health measures are needed to reduce transmission in densely populated areas, particularly in regions with high population density such as cities

``

Download math_model_eco-epidemiology.mw

I'm excited to announce the creation of a new LinkedIn group, Maple Software Community! This group is dedicated to discussions about the use of Maple software and is designed to be a valuable resource for undergraduate and graduate students, researchers, and all Maple enthusiasts.

By joining this community, you'll have the opportunity to:

  • Learn about upcoming events and workshops that can enhance your skills.
  • Stay informed on the latest projects that leverage Maple software.
  • Engage in discussions that explore the many uses of Maple across various fields.
  • Connect with Maple ambassadors and users worldwide who are eager to share their knowledge and experience.

Whether you're a seasoned user or just starting out with Maple, your contributions to this group are welcome and encouraged. Let's build a thriving community together!


Looking forward to seeing you there! 

Maple Software Community

I'm excited to share some valuable resources that I've found incredibly helpful for anyone looking to enhance their Maple skills. Whether you're just starting, studying as a student, or are a seasoned professional, these guides and books offer a wealth of information to aid your learning journey.

Exploring Discrete Mathematics With Maple

These materials are freely available and can be a great addition to your learning resources. They cover a wide range of topics and are designed to help users at all levels improve their Maple proficiency.

You can add your own sources in a comment!

Happy learning and I hope you find these resources as useful as I have!

 Introduction
Maple Coding Expert is a GPT-based AI tool designed to assist with various mathematical tasks using Maple software. It offers step-by-step guidance and detailed explanations for a range of functions, making it a valuable resource for students, educators, and professionals.

 Core Features and Functions

1.Graph Creation:

   - Function Plotting: Users can plot a wide range of mathematical functions. For instance, to plot the function y = x2, the user would input the command `plot(x^2, x = -10..10);` in Maple. The expert helps in setting up the plotting parameters to visualize the function effectively.
   - Advanced Graphing: Beyond simple functions, the expert can guide users through plotting more complex functions and customizing plots with labels, legends, and different styles.

2. Equation Definition and Manipulation:

   - Defining Equations: The tool assists in defining equations for various calculus operations. For example, to differentiate a function, the command might be `diff(f(x), x);`. This helps in accurately modeling the equations necessary for solving real-world problems.
   - Solving Integrals: For integral calculus, users can get assistance in setting up both definite and indefinite integrals. Commands like `int(f(x), x);` are used to perform integration in Maple.

3. Calculus Problem Solving:
   - Differentiation and Integration: The expert provides guidance on solving derivatives and integrals, which are fundamental operations in calculus. It supports both symbolic and numerical methods, allowing users to choose the best approach for their problem.
   - Differential Equations: Users can solve ordinary and partial differential equations using commands like `dsolve({equations}, {variables});`. The expert offers advice on choosing solution methods and interpreting results.

I recently tried using the Maple Coding Expert for solving some calculus problems. It worked well overall and provided detailed solutions, though sometimes it approached the problems in a more complicated way than expected. Despite this, the accuracy and depth of the explanations were impressive and very helpful for understanding the underlying concepts.

 

Maple Coding Expert stands out as a comprehensive tool for anyone involved with Maple software for mathematical computing. It enhances learning, supports professional tasks, and aids in solving complex mathematical problems with ease.

For more information, you can explore the Maple Coding Expert on [GPTs Hunter](https://www.gptshunter.com/gpt-store/MzExMzI2MzYyMzJlNTAxMjM2) and [YesChat.ai](https://www.yeschat.ai).

 

 

Dear Maple Community,

It has been a year since the passing of Stefan Vorkoetter, who started contributing to the Maple project in the 80s and was a long term member of our development team. 

Here are a few recently published articles about Stefan, that I'd like to share with you:

https://mapletransactions.org/index.php/maple/article/view/18269
https://mapletransactions.org/index.php/maple/article/view/18681

we shall not forget

 

 

In our recent project, we're diving deep into understanding the SIR model—a fundamental framework in epidemiology that helps us analyze how diseases spread through populations. The SIR model categorizes individuals into three groups: Susceptible (S), Infected (I), and Recovered (R). By tracking how people move through these categories, we can predict disease dynamics and evaluate interventions.

Key Points of the SIR Model:

  • Susceptible (S): Individuals who can catch the disease.
  • Infected (I): Those currently infected and capable of spreading the disease.
  • Recovered (R): Individuals who have recovered and developed immunity.

Vaccination Impact: One of the critical interventions in disease control is vaccination, which moves individuals directly from the susceptible to the recovered group. This simple action reduces the number of people at risk, thereby lowering the overall spread of the disease.

We're experimenting with a simple model to understand how different vaccination rates can significantly alter the dynamics of an outbreak. By simulating scenarios with varying vaccination coverage, students can observe how herd immunity plays a crucial role in controlling diseases. Our goal is to make these abstract concepts clear and relatable through practical modeling exercises.


 

In this exercise, we are going back to the simple SIR model, without births or deaths, to look at the effect of vaccination. The aim of this activity is to represent vaccination in a very simple way - we are assuming it already happened before we run our model! By changing the initial conditions, we can prepare the population so that it has received a certain coverage of vaccination.

We are starting with the transmission and recovery parameters  b = .4/daysand c = .1/days . To incorporate immunity from vaccination in the model, we assume that a proportion p of the total population starts in the recovered compartment, representing the vaccine coverage and assuming the vaccine is perfectly effective. Again, we assume the epidemic starts with a single infected case introduced into the population.​
We are going to model this scenario for a duration of 2 years, assuming that the vaccine coverage is 50%, and plot the prevalence in each compartment over time.

 

restart
with(plots)

b := .4; c := .1; n := 10^6; p := .5

deS := diff(S(t), t) = -b*S(t)*I0(t); deI := diff(I0(t), t) = b*S(t)*I0(t)-c*I0(t); deR := diff(R(t), t) = c*I0(t)

diff(R(t), t) = .1*I0(t)

(1)

F := dsolve([deS, deI, deR, S(0) = 1-p, I0(0) = 1/n, R(0) = p], [S(t), I0(t), R(t)], numeric, method = rkf45, maxfun = 100000)

odeplot(F, [[t, S(t)], [t, I0(t)], [t, R(t)]], t = 0 .. 730, colour = [blue, red, green], legend = ["S(t)", "I0(t)", "R(t)"], labels = ["Time (days)", "  Proportion\nof Population "], title = "SIR Model with vaccine coverage 50 %", size = [500, 300])

 

F(100)

[t = 100., S(t) = HFloat(0.46146837378273076), I0(t) = HFloat(0.018483974421123688), R(t) = HFloat(0.5200486517961457)]

(2)

eval(S(:-t), F(100))

HFloat(0.46146837378273076)

(3)

Reff := proc (s) options operator, arrow; b*(eval(S(:-t), F(s)))/(c*n) end proc; Reff(100)

HFloat(1.845873495130923e-6)

(4)

plot(Reff, 0 .. 730, size = [500, 300])

 

Increasing the vaccine coverage to 75%

NULL

restart
with(plots)

b := .4; c := .1; n := 10^6; p := .75

deS := diff(S(t), t) = -b*S(t)*I0(t); deI := diff(I0(t), t) = b*S(t)*I0(t)-c*I0(t); deR := diff(R(t), t) = c*I0(t)

diff(R(t), t) = .1*I0(t)

(5)

NULL

F1 := dsolve([deS, deI, deR, S(0) = 1-p, I0(0) = 1/n, R(0) = p], [S(t), I0(t), R(t)], numeric, method = rkf45, maxfun = 100000)

odeplot(F1, [[t, S(t)], [t, I0(t)], [t, R(t)]], t = 0 .. 730, colour = [blue, red, green], legend = ["S(t)", "I0(t)", "R(t)"], labels = ["Time (days)", "  Proportion\nof Population "], title = "SIR Model with vaccine coverage 75%", size = [500, 300])

 

F(1100)

eval(S(:-t), F1(100))

HFloat(0.249990000844159)

(6)

Reff := proc (s) options operator, arrow; b*(eval(S(:-t), F1(s)))/(c*n) end proc; Reff(100)

HFloat(9.99960003376636e-7)

(7)

plot(Reff, 0 .. 730, size = [500, 300])

 

Does everyone in the population need to be vaccinated in order to prevent an epidemic?What do you observe if you model the infection dynamics with different values for p?

No, not everyone in the population needs to be vaccinated in order to prevent an epidemic . In this scenario, if p equals 0.75 or higher, no epidemic occurs - 75 % is the critical vaccination/herd immunity threshold . Remember,, herd immunity describes the phenomenon in which there is sufficient immunity in a population to interrupt transmission . Because of this, not everyone needs to be vaccinated to prevent an outbreak .

What proportion of the population needs to be vaccinated in order to prevent an epidemic if b = .4and c = .2/days? What if b = .6 and "c=0.1 days^(-1)?"

In the context of the SIR model, the critical proportion of the population that needs to be vaccinated in order to prevent an epidemic is often referred to as the "herd immunity threshold" or "critical vaccination coverage."

• 

Scenario 1: b = .4and c = .2/days

``

restart
with(plots)

b := .4; c := .2; n := 10^6; p := .5``

deS := diff(S(t), t) = -b*S(t)*I0(t); deI := diff(I0(t), t) = b*S(t)*I0(t)-c*I0(t); deR := diff(R(t), t) = c*I0(t)

diff(R(t), t) = .2*I0(t)

(8)

F1 := dsolve([deS, deI, deR, S(0) = 1-p, I0(0) = 1/n, R(0) = p], [S(t), I0(t), R(t)], numeric, method = rkf45, maxfun = 100000)

odeplot(F1, [[t, S(t)], [t, I0(t)], [t, R(t)]], t = 0 .. 730, colour = [blue, red, green], legend = ["S(t)", "I0(t)", "R(t)"], labels = ["Time (days)", "  Proportion\nof Population "], title = "SIR Model with vaccine coverage 50 %", size = [500, 300])

 


The required vaccination coverage is around 50% .

• 

Scenario 1: b = .6and c = .1/days

restart
with(plots)

b := .6; c := .1; n := 10^6; p := .83NULL

deS := diff(S(t), t) = -b*S(t)*I0(t); deI := diff(I0(t), t) = b*S(t)*I0(t)-c*I0(t); deR := diff(R(t), t) = c*I0(t)

diff(R(t), t) = .1*I0(t)

(9)

NULL

F1 := dsolve([deS, deI, deR, S(0) = 1-p, I0(0) = 1/n, R(0) = p], [S(t), I0(t), R(t)], numeric, method = rkf45, maxfun = 100000)

odeplot(F1, [[t, S(t)], [t, I0(t)], [t, R(t)]], t = 0 .. 730, colour = [blue, red, green], legend = ["S(t)", "I0(t)", "R(t)"], labels = ["Time (days)", "  Proportion\nof Population "], title = "SIR Model with vaccine coverage 83% ", size = [500, 300])

 

"The required vaccination coverage is around 83 `%` ."


Download SIR_simple_vaccination_example.mw

This Maplesoft guest blog post is from Prof. Dr. Johannes Blümlein from Deutsches Elektronen-Synchrotron (DESY), one of the world’s leading particle accelerator centres used by thousands of researchers from around the world to advance our knowledge of the microcosm. Prof. Dr. Blümlein is a senior researcher in the Theory Group at DESY, where he and his team make significant use of Maple in their investigations of high energy physics, as do other groups working in Quantum Field Theory. In addition, he has been involved in EU programs that give PhD students opportunities to develop their Maple programming skills to support their own research and even expand Maple’s support for theoretical physics.


 

The use of Maple in solving frontier problems in theoretical high energy physics

For several decades, progress in theoretical high energy physics relies on the use of efficient computer-algebra calculations. This applies both to so-called perturbative calculations, but also to non-perturbative computations in lattice field theory. In the former case, large classes of Feynman diagrams are calculated analytically and are represented in terms of classes of special functions. In early approaches started during the 1960s, packages like Reduce [1] and Schoonship [2] were used. In the late 1980s FORM [3] followed and later on more general packages like Maple and Mathematica became more and more important in the solution of these problems. Various of these problems are related to data amounts in computer-algebra of O(Tbyte) and computation times of several CPU years currently, cf. [4].

Initially one has to deal with huge amounts of integrals. An overwhelming part of them is related by Gauss’ divergence theorem up to a much smaller set of the so-called master integrals (MIs). One performs first the reduction to the MIs which are special multiple integrals. No general analytic integration procedures for these integrals exist. There are, however, several specific function spaces, which span these integrals. These are harmonic polylogarithms, generalized harmonic polylogarithms, root-valued iterated integrals and others. For physics problems having solutions in these function spaces codes were designed to compute the corresponding integrals. For generalized harmonic polylogarithms there is a Maple code HyperInt [5] and other codes [6], which have been applied in the solution of several large problems requiring storage of up to 30 Gbyte and running times of several days. In the systematic calculation of special numbers occurring in quantum field theory such as the so-called β-functions and anomalous dimensions to higher loop order, e.g. 7–loop order in Φ4 theory, the Maple package HyperLogProcedures [7] has been designed. Here the largest problems solved require storage of O(1 Tbyte) and run times of up to 8 months. Both these packages are available in Maple only.

A very central method to evaluate master integrals is the method of ordinary differential equations. In the case of first-order differential operators leading up to root-valued iterative integrals their solution is implemented in Maple in [8] taking advantage of the very efficient differential equation solvers provided by Maple. Furthermore, the Maple methods to deal with generating functions as e.g. gfun, has been most useful here. For non-first order factorizing differential equation systems one first would like to factorize the corresponding differential operators [9]. Here the most efficient algorithms are implemented in Maple only. A rather wide class of solutions is related to 2nd order differential equations with more than three singularities. Also here Maple is the only software package which provides to compute the so-called 2F1 solutions, cf. [10], which play a central role in many massive 3-loop calculations

The Maple-package is intensely used also in other branches of particle physics, such as in the computation of next-to-next-to leading order jet cross sections at the Large Hadron Collider (LHC) with the package NNLOJET and double-parton distribution functions. NNLOJET uses Maple extensively to build the numerical code. There are several routines to first build the driver with automatic links to the matrix elements and subtraction terms, generating all of the partonic subprocesses with the correct factors. To build the antenna subtraction terms, a meta-language has been developed that is read by Maple and converted into calls to numerical routines for the momentum mappings, calls to antenna and to routines with experimental cuts and plotting routines, cf. [11].

In lattice gauge calculations there is a wide use of Maple too. An important example concerns the perturbative predictions in the renormalization of different quantities. Within different European training networks, PhD students out of theoretical high energy physics and mathematics took the opportunity to take internships at Maplesoft for several months to work on parts of the Maple package and to improve their programming skills. In some cases also new software solutions could be obtained. Here Maplesoft acted as industrial partner in these academic networks.

References

[1] A.C. Hearn, Applications of Symbol Manipulation in Theoretical Physics, Commun. ACM 14 No. 8, 1971.

[2] M. Veltman, Schoonship (1963), a program for symbolic handling, documentation, 1991, edited by D.N. Williams.

[3] J.A.M. Vermaseren, New features of FORM, math-ph/0010025.

[4] J. Blümlein and C. Schneider, Analytic computing methods for precision calculations in quantum field theory, Int. J. Mod. Phys. A 33 (2018) no.17, 1830015 [arXiv:1809.02889 [hep-ph]].

[5] E. Panzer, Algorithms for the symbolic integration of hyperlogarithms with applications to Feynman integrals, Comput. Phys. Commun. 188 (2015) 148–166 [arXiv:1403.3385 [hep-th]].

[6] J. Ablinger, J. Blümlein, C .Raab, C. Schneider and F. Wissbrock, Calculating Massive 3-loop Graphs for Operator Matrix Elements by the Method of Hyperlogarithms, Nucl. Phys. 885 (2014) 409-447 [arXiv:1403.1137 [hep-ph]].

[7] O. Schnetz, φ4 theory at seven loops, Phys. Rev. D 107 (2023) no.3, 036002 [arXiv: 2212.03663 [hep-th]].

[8] J. Ablinger, J. Blümlein, C. G. Raab and C. Schneider, Iterated Binomial Sums and their Associated Iterated Integrals, J. Math. Phys. 55 (2014) 112301 [arXiv:1407.1822 [hep-th]].

[9] M. van Hoeij, Factorization of Differential Operators with Rational Functions Coefficients, Journal of Symbolic Computation, 24 (1997) 537–561.

[10] J. Ablinger, J. Blümlein, A. De Freitas, M. van Hoeij, E. Imamoglu, C. G. Raab, C. S. Radu and C. Schneider, Iterated Elliptic and Hypergeometric Integrals for Feynman Diagrams, J. Math. Phys. 59 (2018) no.6, 062305 [arXiv:1706.01299 [hep-th]].

[11] A. Gehrmann-De Ridder, T. Gehrmann, E.W.N. Glover, A. Huss and T.A. Morgan, Precise QCD predictions for the production of a Z boson in association with a hadronic jet, Phys. Rev. Lett. 117 (2016) no.2, 022001 [arXiv:1507.02850 [hep-ph]].

To celebrate this day of mathematics, I want to share my favourite equation involving Pi, the Bailey–Borwein–Plouffe (BBP) formula:

This is my favourite for a number of reasons. Firstly, Simon Plouffe and the late Peter Borwein (two of the authors that this formula is named after) are Canadian! While I personally have nothing to do with this formula, the fact that fellow Canadians contributed to such an elegant equation is something that I like to brag about.

Secondly, I find it fascinating how Plouffe first discovered this formula using a computer program. It can often be debated whether mathematics is discovered or invented, but there’s no doubt here since Plouffe found this formula by doing an extensive search with the PSLQ integer relation algorithm (interfaced with Maple). This is an example of how, with ingenuity and creativity, one can effectively use algorithms and programs as powerful tools to obtain mathematical results.

And finally (most importantly), with some clever rearranging, it can be used to compute arbitrary digits of Pi!

Digit 2024 is 8
Digit 31415 is 5
Digit 123456 is 4
Digit 314159 is also 4
Digit 355556 is… F?

That last digit might look strange… and that’s because they’re all in hexadecimal (base-16, where A-F represent 10-15). As it turns out, this type of formula only exists for Pi in bases that are powers of 2. Nevertheless, with the help of a Maple script and an implementation of the BBP formula by Carl Love, you can check out this Learn document to calculate some arbitrary digits of Pi in base-16 and learn a little bit about how it works.

After further developments, this formula led to project PiHex, a combined effort to calculate as many digits of Pi in binary as possible; it turns out that the quadrillionth bit of Pi is zero! This also led to a class of BBP-type formulas that can calculate the digits of other constants like (log2)*(π^2) and (log2)^5.

Part of what makes this formula so interesting is human curiosity: it’s fun to know these random digits. Another part is what makes mathematics so beautiful: you never know what discoveries this might lead to in the future. Now if you’ll excuse me, I have a slice of lemon meringue pie with my name on it 😋

 

References
BBP Formula (Wikipedia)
A Compendium of BBP-Type Formulas
The BBP Algorithm for Pi

In the most recent issue of Maple Transactions, I published (with David Jeffrey, and with a student named Johan Joby) a paper that used Jupyter Notebook with a Maple kernel as the main vehicle.  Have a look, and let me know what you think.

Two-cycles in the infinite exponential tower

1 2 3 4 5 6 7 Last Page 1 of 26