Applications, Examples and Libraries

Share your work here


This post is inspired by minhthien2016's question.

The problem, denoted 2/N/1, for reasons that will appear clearly further on, is to pack N disks into the unit square in such a way that the sum of their radii is maximum.

I replied this problem using Optimization-NLPSolve for N from 1 (obvious solution) to 16, which motivated a few questions, in particular:

  • @Carl Love: "Can we confirm that the maxima are global (NLPSolve tends to return local optima)?
    Using NLPSolve indeed does not guarantee that the solution found is the (a?) global maximum. In fact packing problems are generaly tackled by using evolutionnary algorithms, greedy algorithms, or specific heuristic strategies.
    Nevertheless, running NLPSolve a large number of times from different initial points may provide several different optima whose the largest one can be reasonably considered as the global maximum you are looking for.
    Of course this may come to a large price in term of computational time.

     
  • @acer: "How tight are [some constraints], at different N? Are they always equality?"
    The fact some inequality constraints type always end to equality constraints (which means that in an optimal packing each disk touches at least one other annd, possibly the boundary of the box) seems hard to prove mathematically, but I gave here a sketch of informal proof.



I found 2/N/1 funny enough to spend some time digging into it and looking to some generalizations I will refer to as D/N/M:  How to pack N D-hypersheres into the unit D-hypercube such that the sum of the M-th power of their radii is maximum?
For the sake of simplicity I will say ball instead of disk/sphere/hypersphere and box instead of square/cube/hypercube.

The first point is that problems D/N/1 do not have a unique solution as soon as N > 1 , indeed any solution can be transformed into another one using symmetries with respect to medians and diagonals of the box. Hereafter I use this convention:

Two solutions and s' are fundamental solutions if:

  1. the ordered lists of radii and s'  contain are identical but there is no composition of symmetries from to s',
  2. or, the ordered lists of radii and s'  contain are not identical.
     

It is easy to prove that 2/2/1 and 3/2/1, and likely D/2/1, have an infinity of fundamental solutions: see directory FOCUS__2|2|1_and_3|2|1 in the attached zip file..
At the same time 2/N/2, 3/N/3, and likely D/N/D, have only one fundamental solution (see directory FOCUS__2|N|2 for more details and a simple way to characterize these solutions

 (Indeed the strategy ito find the solution of D/N/D  in placing the biggest possible ball in the largest void D/N-1/D contains. Unfortunately this characterization is not algorithmically constructive in the sense that findind this biggest void is a very complex geometrical and combinatorial problem.
 it requires finding the largest void  in a pack of balls)


Let Md, 1(N)  the maximum value of the sum of balls radii for problem d/N/1.
The first question I asked myself is: How does Md, 1(N) grows with N?

 

(Md, 1(N) is obviously a strictly increasing function of N: indeed the solution of problem d/N/1 contains several voids where a ball of strictly positive radius can be placed, then  Md+1, 1(N) > Md, 1(N) )


The answer seems amazing as intensive numerical computations suggest that
                                      

See D|N|M__Growth_law.mw in the attached sip file.
This formula fits very well the set of points  { [n, Sd, 1(n) , n=1..48) } for d=2..6.
I have the feeling that this conjecture might be proven (rejected?) by rigourous mathematical arguments.


Fundamental solutions raise several open problems:

  • Are D/2/1 problems the only one with more than one fundamental solutions?

    The truth is that I have not been capable to find any other example (which does not mean they do not exist).
    A quite strange thing is the behaviour of NLPSolve: as all the solutions of D/2/1 are equally likely, the histogram of the solutions provided by a large number of NLPSolve runs from different initial points is far from being uniform.
    F
    or more detail refer ro directory FOCUS__2|2|1_and_3|2|1
     in the attached zip file
    I do not understand where this bias comes from: is it due to the implementation of SQP in NLPSolve, or to SQP itself?

     
  • For some couples (D, N) the solution of D/N/1 is made of balls of same radius.
    For N from 1 to 48 this is (numerically)
     the case for 2/1/1 and2/2/1, but the three dimensional case is reacher as it contains  3/1/13/2/1,  3/3/1,  3/4/1 and 3/14/1 (this latter being quite surprising).
    Is there only a finite number of values 
    N such that D/N/1 is made of balls with identical radii?
    If it is so, is this number increasing with
     D?
    It is worth noting that those values of
    N mean that the solution of problems D/N/1 are identical to those of a more classic packing problem: "What is the largest radius N identical balls packed in a unit bow may have?".
    For an exhaustive survey of this latter problem see
    Packomania.

     
  • A related question is "How does the number of different radii evolves as N increases dor given values of D and M?
    Displays of 2D and 3D packings show that the set of radii has significantly less elements than
    N... at least for values of N not too large. So might we expect that solution of, let us say, 2/100/1 can contain 100 balls of 10 different radii, or it is more reasonable to expect it contains 100 balls of 100 different radii?

     
  • At the opposite numerical investigations of  2/N/1 and  3/N/1 suggest that the number of different radii a fundamental solution contains increases with N (more a trend than a continuous growth).
    So, is it true that very large values of N correspond to solutions where the number of different radii is also very large?

    Or could it be that the growth of the number of different radii I observed is simply the consequence of partially converged results?
     
  • Numerical investigations show that for a given dimension d and a given number of balls n,  solutions of d/n/1 and d/n/M (1 < M < d) problems are rather often the same. Is this a rule with a few exceptions or a false impression due to the fact that I did not pushed the simulations to values of n large enough to draw a more solid picture)?


It is likely that some of the questions above could be adressed by using a more powerful machine than mine.


All the codes and results are gathered in  a zip file you can download from OneDrive Google  (link at the end of this post, 262 Mb, 570 Mb when unzipped, 1119 files).
Install this zip file in the directory of your choice and unzip it to get a directory named
PACKING
Within it:

  • README.mw contains a description of the different codes and directories
  • Repository.rtf must contain a string repesenting the absolute path of directory PACKING


Follow this link OneDrive Google


 

An attractor is called strange if it has a fractal structure, that is if it has a non-integer Hausdorff dimension. This is often the case when the dynamics on it are chaotic, but strange nonchaotic attractors also exist.  If a strange attractor is chaotic, exhibiting sensitive dependence on initial conditions, then any two arbitrarily close alternative initial points on the  attractor, after any of various numbers of iterations, will lead to  points that are arbitrarily far apart (subject to the confines of the attractor), and after any of various other numbers of iterations will  lead to points that are arbitrarily close together. Thus a dynamic system with a chaotic attractor is locally unstable yet globally stable: once some sequences have entered the attractor, nearby points diverge  from one another but never depart from the attractor.


The term strange attractor was coined by David Ruelle and Floris Takens to describe the attractor resulting from a series of bifurcations of a system describing fluid flow. Strange attractors are often differentiable in a few directions, but some are like a Cantor dust, and therefore not differentiable. Strange attractors may also be found  in the presence of noise, where they may be shown to support invariant  random probability measures of Sinai–Ruelle–Bowen type.


Examples of strange attractors include the  Rössler attractor, and Lorenz attractor.

 

 

THOMAS``with(plots); b := .20; sys := diff(x(t), t) = sin(y(t))-b*x(t), diff(y(t), t) = sin(z(t))-b*y(t), diff(z(t), t) = sin(x(t))-b*z(t); sol := dsolve({sys, x(0) = 1.1, y(0) = 1.1, z(0) = -0.1e-1}, {x(t), y(t), z(t)}, numeric); odeplot(sol, [x(t), y(t), z(t)], t = 0 .. 600, axes = boxed, numpoints = 50000, labels = [x, y, z], title = "Thomas Attractor")

 

 

 

Dabras``

with(plots); a := 3.00; b := 2.7; c := 1.7; d := 2.00; e := 9.00; sys := diff(x(t), t) = y(t)-a*x(t)+b*y(t)*z(t), diff(y(t), t) = c*y(t)-x(t)*z(t)+z(t), diff(z(t), t) = d*x(t)*y(t)-e*z(t); sol := dsolve({sys, x(0) = 1.1, y(0) = 2.1, z(0) = -2.00}, {x(t), y(t), z(t)}, numeric); odeplot(sol, [x(t), y(t), z(t)], t = 0 .. 100, axes = boxed, numpoints = 35000, labels = [x, y, z], title = "Dabras Attractor")

 

Halvorsen

NULLwith(plots); a := 1.89; sys := diff(x(t), t) = -a*x(t)-4*y(t)-4*z(t)-y(t)^2, diff(y(t), t) = -a*y(t)-4*z(t)-4*x(t)-z(t)^2, diff(z(t), t) = -a*z(t)-4*x(t)-4*y(t)-x(t)^2; sol := dsolve({sys, x(0) = -1.48, y(0) = -1.51, z(0) = 2.04}, {x(t), y(t), z(t)}, numeric, maxfun = 300000); odeplot(sol, [x(t), y(t), z(t)], t = 0 .. 600, axes = boxed, numpoints = 35000, labels = [x, y, z], title = "Halvorsen Attractor")

 

Chen

 

 

with(plots); alpha := 5.00; beta := -10.00; delta := -.38; sys := diff(x(t), t) = alpha*x(t)-y(t)*z(t), diff(y(t), t) = beta*y(t)+x(t)*z(t), diff(z(t), t) = delta*z(t)+(1/3)*x(t)*y(t); sol := dsolve({sys, x(0) = -7.00, y(0) = -5.00, z(0) = -10.00}, {x(t), y(t), z(t)}, numeric); odeplot(sol, [x(t), y(t), z(t)], t = 0 .. 100, axes = boxed, numpoints = 35000, labels = [x, y, z], title = "Chen Attractor")

 

References

1. 

https://www.dynamicmath.xyz/strange-attractors/

2. 

https://en.wikipedia.org/wiki/Attractor#Strange_attractor

``


 

Download Attractors.mw


 

This project discusses predator-prey system, particularly the Lotka-Volterra equations,which model the interaction between two sprecies: prey and predators. Let's solve the Lotka-Volterra equations numerically and visualize the results.

NULL

NULL

alpha := 1.0; beta := .1; g := 1.5; delta := 0.75e-1; ode1 := diff(x(t), t) = alpha*x(t)-beta*x(t)*y(t); ode2 := diff(y(t), t) = delta*x(t)*y(t)-g*y(t); eq1 := -beta*x*y+alpha*x = 0; eq2 := delta*x*y-g*y = 0; equilibria := solve({eq1, eq2}, {x, y}); print("Equilibrium Points: ", equilibria); initial_conditions := x(0) = 40, y(0) = 9; sol := dsolve({ode1, ode2, initial_conditions}, {x(t), y(t)}, numeric); eq_points := [seq([rhs(eq[1]), rhs(eq[2])], `in`(eq, equilibria))]

[[0., 0.], [20., 10.]]

plots[odeplot](sol, [[t, x(t)], [t, y(t)]], t = 0 .. 100, legend = ["Rabbits", "Wolves"], title = "Prey-Predator Dynamics", labels = ["Time", "Population"])

NULL

NULL

NULL

sol_plot := plots:-odeplot(sol, [[x(t), y(t)]], 0 .. 100, color = "blue"); equilibrium_plot := plots:-pointplot(eq_points, color = "red", symbol = solidcircle, symbolsize = 15); plots:-display([sol_plot, equilibrium_plot], title = "Phase Portrait with Equilibrium Points", labels = ["Rabbits", "Wolves"])

Now, we need to handle a modified version of the Lotka-Volterra equations. These modified equations incorporate logistic growth fot the prey population.

 

 

restart

alpha := 1.0; beta := .1; g := 1.5; delta := 0.75e-1; k := 100; ode1 := diff(x(t), t) = alpha*x(t)*(1-x(t)/k)-beta*x(t)*y(t); ode2 := diff(y(t), t) = delta*x(t)*y(t)-g*y(t); eq1 := alpha*x*(1-x/k)-beta*x*y = 0; eq2 := delta*x*y-g*y = 0; equilibria := solve({eq1, eq2}, {x, y}); print("Equilibrium Points: ", equilibria); initial_conditions := x(0) = 40, y(0) = 9; sol := dsolve({ode1, ode2, initial_conditions}, {x(t), y(t)}, numeric); eq_points := [seq([rhs(eq[1]), rhs(eq[2])], `in`(eq, equilibria))]

[[0., 0.], [100., 0.], [20., 8.]]

plots[odeplot](sol, [[t, x(t)], [t, y(t)]], t = 0 .. 100, legend = ["Rabbits", "Wolves"], title = "Prey-Predator Dynamics", labels = ["Time", "Population"])

NULL

plots:-odeplot(sol, [[x(t), y(t)]], 0 .. 50, color = "blue"); equilibrium_plot := plots:-pointplot(eq_points, color = "red", symbol = solidcircle, symbolsize = 15); plots:-display([plots:-odeplot(sol, [[x(t), y(t)]], 0 .. 50, color = "blue"), equilibrium_plot], title = "Phase Portrait with Equilibrium Points", labels = ["Rabbits", "Wolves"])

NULL


 

Download predator_prey2.mw

We are pleased to announce that the registration for the Maple Conference 2024 is now open.

Like the last few years, this year’s conference will be a free virtual event. Please visit the conference page for more information on how to register.

This year we are offering a number of new sessions, including more product training options and an Audience Choice session.
You can find an overview of the program on the Sessions page. Those who register before September 10th, 2024 will have a chance to vote for the topics they want to learn more about during the Audience Choice session.

We hope to see you there!

This is a reminder that presentation applications for the Maple Conference are due July 17, 2024.

The conference is a a free virtual event and will be held on October 24 and 25, 2024.

We are inviting submissions of presentation proposals on a range of topics related to Maple, including Maple in education, algorithms and software, and applications. We also encourage submission of proposals related to Maple Learn. You can find more information about the themes of the conference and how to submit a presentation proposal at the Call for Participation page.

I encourage all of you here in the Maple Primes community to consider joining us for this event, whether as a presenter or an attendee!

Kaska Kowalska
Contributed Program Co-Chair

 

The Proceedings of the Maple Conference 2023 is now out, at

mapletransactions.org

The presentations these are based on (and more) can be found at https://www.maplesoft.com/mapleconference/2023/full-program.aspx#schedule .

There are several math research papers using Maple, an application paper by an undergraduate student, an engineering application paper, and an interesting geometry teaching paper.

Please have a look, and don't forget to register for the Maple Conference 2024.

Consider the equation  (2^x)*(27^(1/x)) = 24  for which we need to find the exact values ​​of its real roots. This is not difficult to solve by hand if you first take the logarithm of this equation to any base, after which the problem is reduced to solving a quadratic equation. But the  solve  command fails to solve this equation and returns the result in RootOf form. The problem is solved if we first ask Maple to take the logarithm of the equation. I wonder if the latest versions of Maple also do not directly address the problem?

restart;
Eq:=2^x*27^(1/x)=24:
solve(Eq, x, explicit);

map(ln, Eq); # Taking the logarithm of the equation
solve(%, x);
simplify({%}); # The final result

                  

 

We've just launched Maple Flow 2024!

You're in the driving seat with Maple Flow - each new feature has a straight-line connection to a user-driven demand to work faster and more efficiently.

Head on over here for a rundown of everything that's new, but I thought I'd share my personal highlights here.

If your result contains a large vector or matrix, you can now scroll to see more data. You can also change the size of the matrix to view more or fewer rows and columns.

You can resize rows and columns if they're too large or small, and selectively enable row and column headers.

If the vector or matrix in your result contains a unit, you can now rescale units with the Context Panel (for the entire matrix) or inline (for individual entries).

A few releases ago, we introduced the Variables palette to help you keep track of all the user-defined parameters at point of the grid cursor.

You can now insert variables into the worksheet from the Variables palette. Just double-click on the appropriate name.

Maple Flow already features command completion - just type the first few letters of a command, and a list of potential completions appears. Just pick the completion you need with a quick tap of the Tab key.

We've supercharged this feature to give potential arguments for many popular functions. Type a function name followed by an opening bracket, and a list appears.

In case you've missed it, the argument completion list also features (when they make sense) user-defined variables.

You can now link to different parts of the same worksheet. This can be used to create a table of contents that lets you jump to different parts of larger worksheets.

This page lists everything that's new in the current release, and all the prior releases. You might notice that we have three releases a year, each featuring many user-requested items. Let me know what you want to see next - you might not have to wait that long!

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


 

When a derivative can be written as a function of the independent variable only for example

y'=f(x)

y''=f(x)

y'''(x)=f(x)

etc.

 

We call that a directly integrable equation.

 

Example 1:

 

Find the general solution for the following directly integrable equation

diff(y, x) = 6*x^2+4*y(1) and 6*x^2+4*y(1) = 0

That means

int(6*x^2+4, x)

y = 2*x^3+c+4*x", where" c is an arbitary solution

``

 

 
equation1 := diff(y(x), x) = 6*x^2+4

diff(y(x), x) = 6*x^2+4

(1)

NULL

NULL

sol1 := dsolve(equation1, y(x))

y(x) = 2*x^3+c__1+4*x

(2)

And  if we have the initial condition
y(1) = 0
particular_sol1 := dsolve({equation1, y(1) = 0}, y(x))

y(x) = 2*x^3+4*x-6

(3)

"(->)"

 

 

 

 

 

"Example 2:"NULL

NULL

"  Find the particular solution for the following equation with condition"

 

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

y(1)=3

So we will need to get the y' by itself

int(-1/x^2, x)

so,

y = 1/x+c , where c is an arbitary constant

And this is our general solution. Now we plug in the initial condition when x = 1, y = 3.

 

That means c = 2.

 

Thus, the particular solution is

 

y = 1/x+2``

eq := x^2*(diff(y(x), x)) = -1

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

(4)

NULL

NULL

sol := dsolve(eq, y(x))

y(x) = 1/x+c__1

(5)

NULL

particular_sol := dsolve({eq, y(1) = 3}, y(x))

y(x) = 1/x+2

(6)

NULL

NULL

plot(1/x+2, x = -20 .. 20, color = "Red", axes = normal, legend = [typeset(1/x+2)])

 

NULL

NULL

NULL

NULL

" Example 3:"

 

" Find the particular solution for the following equation with condition"

 

diff(y, t, t) = cost, (D(y))(0) = 0, y(0) = 1

eq1 := diff(y(t), t, t) = cos(t)

diff(diff(y(t), t), t) = cos(t)

(7)

particular_sol := dsolve({eq1, y(0) = 1, (D(y))(0) = 0}, y(t))

y(t) = -cos(t)+2

(8)

"(->)"

 

 

 

 

NULL


 

Download integral.mw

 

Hello,

Attached I am sending several procedures for curves in 3D space. They were written without using Maple's built-in DiffGeo procedures and functions. As an example of their use, I made several animations - Maple worksheets are attached. I hope that maybe they will be useful to someone.

Regards.

ClsDGproc-Curves.zip

restart;

torus:= (x^2+y^2+z^2 + R^2-r^2)^2 - 4*R^2*(x^2+y^2):

torus1:=eval(torus,[r=1,R=4]);

(x^2+y^2+z^2+15)^2-64*x^2-64*y^2

(1)

plots:-implicitplot3d(torus1=0, x=-6..6,y=-6..6,z=-6..6, numpoints=5000, scaling=constrained, view=-2..2);

 

sol1:=solve(torus1<0, [x,y,z]); # It should be easy for Maple

[[x = 0, y < -3, -5 < y, z < (-y^2-8*y-15)^(1/2), -(-y^2-8*y-15)^(1/2) < z], [x = 0, y < 5, 3 < y, z < (-y^2+8*y-15)^(1/2), -(-y^2+8*y-15)^(1/2) < z]]

(2)

eval(torus1, [x=3,y=2,z=0]); # in the interior

-48

(3)

eval(sol1, [x=3,y=2,z=0]);   # ???

[[3 = 0, 2 < -3, -5 < 2, 0 < (-35)^(1/2), -(-35)^(1/2) < 0], [3 = 0, 2 < 5, 3 < 2, 0 < (-3)^(1/2), -(-3)^(1/2) < 0]]

(4)

 

 

Download Bug-solve_ineqs.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]].

     Happy Easter to all those who celebrate! One common tradition this time of year is decorating Easter eggs. So, we’ve decided to take this opportunity to create some egg-related math content in Maple Learn. This year, a blog post by Tony Finch inspired us to create a walkthrough exploring the four-point egg. The four-point egg is a method to construct an egg-shaped graph using just a compass and a ruler, or in this case, Maple Learn. Here's the final product: 

     The Maple Learn document, found here, walks through the steps. In general, each part of the egg is an arc corresponding to part of a circle centred around one of the points generated in this construction. 

     For instance, starting with the unit circle and the three red points in the image below, the blue circle is centred at the bottom point such that it intersects with the top of the unit circle, at (0,1). The perpendicular lines were constructed using the three red points, such that they intersect at the bottom point and pass through opposite side points, either (-1,0) or (1,0). Then, the base of the egg is constructed by tracing an arc along the bottom of the blue circle, between the perpendicular lines, shown in red below.

 

     Check out the rest of the steps in the Maple Learn Document. Also, be sure to check out other egg-related Maple Learn documents including John May’s Egg Formulas, illustrating other ways to represent egg-shaped curves with mathematics, and Paige Stone’s Easter Egg Art, to design your own Easter egg in Maple Learn. So, if you’ve had your fill of chocolate eggs, consider exploring some egg-related geometry - Happy Easter!  

Here are the Maple sources for the paper "Maximum Gap among Integers Having a Common Divisor with an Odd semi-prime" published on Journal of Advances in Mathematics and Computer Science 39 (10):51-61, at :https://doi.org/10.9734/jamcs/2024/v39i101934.

 

gap := proc(a, b) return abs(a - b) - 1; end proc

HostsNdivisors := proc(N)

local i, j, g, d, L, s, t, m, p, q, P, Q, np, nq;

m := floor(1/2*N - 1/2);

L := evalf(sqrt(N));

P := Array();

Q := Array();

s := 1; t := 1;

for i from 3 to m do

   d := gcd(i, N);

    if 1 < d and d < L then P(s) := i; s++;

    elif L < d then Q(t) := i; t++; end if;

end do;

  np := s - 1;

  nq := t - 1;

 for i to np do printf("%3d,", P(i)); end do;

  printf("\n");

  for i to nq do printf("%3d,", Q(i)); end do;

  printf("\n gaps: \n");

  for i to np do

     for j to nq do

      p := P(i); q := Q(j);

      g := gap(p, q);

      printf("%4d,", g);

  end do;

    printf("\n");

end do;

end proc

 

HostOfpq := proc(p, q)

local alpha, s, t, g, r, S, T, i, j;

   S := 1/2*q - 1/2;

   T := 1/2*p - 1/2;

   alpha := floor(q/p);

    r := q - alpha*p;

   for s to S do

     for t to T do

       g := abs((t*alpha - s)*p + t*r) - 1;

        printf("%4d,", g);

      end do;

     printf("\n");

 end do;

end proc

MapleSource.pdf

MapleSource.mw

 

3 4 5 6 7 8 9 Last Page 5 of 76