Hi,

I need to run the following procedure a couple of million times. Although it works, Maple sometimes chokes for no apparent reason (if there is a reason, please let me know). I was wondering whether an expert could help me tweak the procedure (or possibly rewrite it) to achieve the best possible performance. I am planning to use Grid:-Map or, if possible, Threads:-Map.

 

generateNonlinearModelsPlus := proc(model::list,fullmodel::list,vars::list:=[x,y,z])
description "This function generates a list of all models with one more monomial from the full model":
local tab::table(),n:=nops(model),i,j,k:=1,ans,terms,aaa,allmoncoefThreads:
# local procedure
allmoncoefThreads := proc(f::list,vars::list)
description "This function finds the monomials multipled by their coefficients for each expression (equation) of a list.":
local n:=numelems(f),i,mon:=[seq](0,i=1..n),M,cc:=[seq](0,i=1..n),ans:
for i from 1 to n do
  cc[i]:=[coeffs](expand(f[i]),vars, 'M'):
  mon[i]:=[M]:
end do:
ans:=[seq](zip((ww,vv)->ww*vv,cc[i],mon[i]),i=1..n):
return(ans)
end proc:
# main part
ans:=zip((w,v)->expand(simplify(v-w)),model,fullmodel): # Find the monomials that are not in model
terms:=allmoncoefThreads(ans,vars): # Separate the monomials
#
for i from 1 to n do
   aaa:=model:
   for j from 1 to nops(terms[i]) do
       aaa[i]:=model[i]+terms[i,j]:
       tab[k]:=aaa:
       k:=k+1:
   end do:
end do:
tab:=convert(tab,list):
return(tab):
end proc:

Here is an example of how I run it: 

model:=[y*alpha[1, 2], z*alpha[2, 3], x^3*alpha[3, 10] + x*alpha[3, 1] + alpha[3, 0]]:

fullmodel:=[x^3*alpha[1, 10] + x^2*y*alpha[1, 11] + x^2*z*alpha[1, 12] + x*y^2*alpha[1, 13] + x*y*z*alpha[1, 14] + x*z^2*alpha[1, 15] + y^3*alpha[1, 16] + y^2*z*alpha[1, 17] + y*z^2*alpha[1, 18] + z^3*alpha[1, 19] + x^2*alpha[1, 4] + x*y*alpha[1, 5] + x*z*alpha[1, 6] + y^2*alpha[1, 7] + y*z*alpha[1, 8] + z^2*alpha[1, 9] + x*alpha[1, 1] + y*alpha[1, 2] + z*alpha[1, 3] + alpha[1, 0], x^3*alpha[2, 10] + x^2*y*alpha[2, 11] + x^2*z*alpha[2, 12] + x*y^2*alpha[2, 13] + x*y*z*alpha[2, 14] + x*z^2*alpha[2, 15] + y^3*alpha[2, 16] + y^2*z*alpha[2, 17] + y*z^2*alpha[2, 18] + z^3*alpha[2, 19] + x^2*alpha[2, 4] + x*y*alpha[2, 5] + x*z*alpha[2, 6] + y^2*alpha[2, 7] + y*z*alpha[2, 8] + z^2*alpha[2, 9] + x*alpha[2, 1] + y*alpha[2, 2] + z*alpha[2, 3] + alpha[2, 0], x^3*alpha[3, 10] + x^2*y*alpha[3, 11] + x^2*z*alpha[3, 12] + x*y^2*alpha[3, 13] + x*y*z*alpha[3, 14] + x*z^2*alpha[3, 15] + y^3*alpha[3, 16] + y^2*z*alpha[3, 17] + y*z^2*alpha[3, 18] + z^3*alpha[3, 19] + x^2*alpha[3, 4] + x*y*alpha[3, 5] + x*z*alpha[3, 6] + y^2*alpha[3, 7] + y*z*alpha[3, 8] + z^2*alpha[3, 9] + x*alpha[3, 1] + y*alpha[3, 2] + z*alpha[3, 3] + alpha[3, 0]]:

vars:=[x,y,z]:

ans:=generateNonlinearModelsPlus(model,fullmodel,vars)

Many thanks.

For quite some time, I have wanted to solve the system attached in "test" using Maple. The smallest solution in natural numbers x, y, and z test.mw

>  (1) >  (2) (3) (4)

Download test.mw

is known, and all these numbers are less than 4 × 10¹². Is this possible in Maple?

(x=1633780814400; y=252782198228; z=3474741058973)

Currently I have Maple versions 2023,2025, and 2026 installed on Windows 11. Today I installed a workbook package containing a module that I just completed using the PackageTools installer in Maple 2026.. To my surprise, I found that a package installed from Maple 2026 was also available in Maple 2023 and, conversely, a package installed in Maple 2023 was automatically available in Maple 2026. i noticed that, with the exception of the Maple Customer Support Updates, the toolbox directory is no longer broken down by versions. I also noticed that the directory containing the module installed by Maple 2026 was named by the workbook instead of the module name (ie. hopfwords.maple). As I recall, the toolboxes used to be version dependent. 

The question is to what extent can one assume that a package created in Maple 2026 will be compatible with at least the more recent versions of Maple, I am also wondering why the directory name is now the workbook name instead of the module name. 

restart;

with(plots): with(LinearAlgebra):

 

# TFSB Coefficients (symbolic in u)

beta0 := u -> (sin(u)*u^3 - 12*u^2 - 24*cos(u) + 24)/(12*(sin(u)*u + 2*cos(u) - 2)*u^2):

beta1 := u -> (5*sin(u)*u^3 + 12*cos(u)*u^2 + 24*cos(u) - 24)/(6*(sin(u)*u + 2*cos(u) - 2)*u^2):

beta2 := u -> beta0(u):

rho0 := u -> ((-u^2-12)*cos(u) - 5*u^2 + 12)/(12*(sin(u)*u + 2*cos(u) - 2)*u^2):

rho1 := u -> (-7*cos(u)*u^3 + 27*sin(u)*u^2 + 120*sin(u) - 120*u)/(60*u^2*(cos(u)*u + 2*u - 3*sin(u))):

rho2 := u -> -rho0(u):

 

# Secondary coefficients (simplified versions)

beta00 := u -> 13/42 - 9*u^2/7840:

beta10 := u -> 1/6 + u^2/720:

beta20 := u -> 1/42 - 17*u^2/70560:

beta01 := u -> 187/1680 + 611*u^2/705600:

beta11 := u -> 11/30 - 29*u^2/25200:

beta21 := u -> 37/1680 + 67*u^2/235200:

beta02 := u -> 11/70 + 491*u^2/352800:

beta12 := u -> 9/10 - 31*u^2/8400:

beta22 := u -> 31/70 + 811*u^2/352800:

 

rho01 := u -> 2/105 + 407*u^2/1058400:

rho11 := u -> -19/210 + 41*u^2/105840:

rho21 := u -> -1/168 - 101*u^2/529200:

rho02 := u -> 53/1680 + 1633*u^2/2116800:

rho12 := u -> 8/105 - 4*u^2/6615:

rho22 := u -> -101/1680 - 2273*u^2/2116800:

 

# Problem definition

omega := 1:

epsilon := 3*Pi/2:

phi := x -> 3*sin(x) - 5*cos(x):  # history function

 

f := (x, v, vp, vd) -> -v - vd + 3*cos(x) + 5*sin(x):

g := proc(x, v, vp, vd, vdp)

    local fx, fv, fvp, fvd;

    fx := -3*sin(x) + 5*cos(x);

    fv := -1;

    fvp := 0;

    fvd := -1;

    return fx + fv*vp + fvp*0 + fvd*vdp;

end proc:

 

# Initial conditions

a := 0: b := 10:

v0 := -5: vp0 := 3:

 

# Variable step-size parameters

tol := 1e-10:

h_min := 0.01:

h_max := 0.5:

h_init := Pi/8:

 

# Store results

X := [a]: V := [v0]: Vp := [vp0]:

h_curr := h_init:

x_curr := a:

v_curr := v0:

vp_curr := vp0:

 

# For history: need v at x-epsilon

get_v_delayed := proc(xx)

    if xx < a then return phi(xx);

    else

        # Interpolate from stored solution

        idx := 1;

        while idx < nops(X) and X[idx] < xx do idx := idx+1; end do;

        if idx = 1 then return phi(xx);

        elif X[idx] = xx then return V[idx];

        else

            # Linear interpolation

            return V[idx-1] + (V[idx]-V[idx-1])*(xx-X[idx-1])/(X[idx]-X[idx-1]);

        end if;

    end if;

end proc:

 

# Newton solver for block

solve_block := proc(x0, v0, vp0, h, omega)

    local u, bet0, bet1, bet2, rho0, rho1, rho2, bet00, bet10, bet20, bet01, bet11, bet21, bet02, bet12, bet22,

          rho01, rho11, rho21, rho02, rho12, rho22, F, J, V0, V1, V2, Vp0, Vp1, Vp2, tolN, iter, dv, dV;

   

    u := omega*h;

    bet0 := beta0(u); bet1 := beta1(u); bet2 := beta2(u);

    rho0 := rho0(u); rho1 := rho1(u); rho2 := rho2(u);

    bet00 := beta00(u); bet10 := beta10(u); bet20 := beta20(u);

    bet01 := beta01(u); bet11 := beta11(u); bet21 := beta21(u);

    bet02 := beta02(u); bet12 := beta12(u); bet22 := beta22(u);

    rho01 := rho01(u); rho11 := rho11(u); rho21 := rho21(u);

    rho02 := rho02(u); rho12 := rho12(u); rho22 := rho22(u);

   

    # Initial guesses

    V1 := v0 + h*vp0;

    V2 := v0 + 2*h*vp0;

    Vp1 := vp0;

    Vp2 := vp0;

   

    tolN := 1e-12;

    for iter from 1 to 10 do

        # Compute delayed values

        vd0 := get_v_delayed(x0 - epsilon);

        vd1 := get_v_delayed(x0 + h - epsilon);

        vd2 := get_v_delayed(x0 + 2*h - epsilon);

        vdp0 := (get_v_delayed(x0 - epsilon + 1e-8) - vd0)/1e-8;

        vdp1 := (get_v_delayed(x0 + h - epsilon + 1e-8) - vd1)/1e-8;

        vdp2 := (get_v_delayed(x0 + 2*h - epsilon + 1e-8) - vd2)/1e-8;

       

        # Compute gamma and g

        gam0 := f(x0, v0, vp0, vd0);

        gam1 := f(x0+h, V1, Vp1, vd1);

        gam2 := f(x0+2*h, V2, Vp2, vd2);

        g0 := g(x0, v0, vp0, vd0, vdp0);

        g1 := g(x0+h, V1, Vp1, vd1, vdp1);

        g2 := g(x0+2*h, V2, Vp2, vd2, vdp2);

       

        # Residuals

        F1 := h*vp0 - (V1 - v0 + h^2*(bet00*gam0 + bet10*gam1 + bet20*gam2)

              + h^3*(rho01*g0 + rho11*g1 + rho21*g2));

        F2 := h*Vp1 - (V1 - v0 + h^2*(bet01*gam0 + bet11*gam1 + bet21*gam2)

              + h^3*(rho01*g0 + rho11*g1 + rho21*g2));

        F3 := h*Vp2 - (V1 - v0 + h^2*(bet02*gam0 + bet12*gam1 + bet22*gam2)

              + h^3*(rho02*g0 + rho12*g1 + rho22*g2));

        F4 := V2 - (2*V1 - v0 + h^2*(bet0*gam0 + bet1*gam1 + bet2*gam2)

              + h^3*(rho0*g0 + rho1*g1 + rho2*g2));

       

        F := Vector([F1, F2, F3, F4]);

        if LinearAlgebra:-Norm(F) < tolN then break; end if;

       

        # Approximate Jacobian (finite differences)

        J := Matrix(4,4);

        delta := 1e-6;

        for j from 1 to 4 do

            V_pert := Vector([V1, V2, Vp1, Vp2]);

            V_pert[j] := V_pert[j] + delta;

            V1p := V_pert[1]; V2p := V_pert[2]; Vp1p := V_pert[3]; Vp2p := V_pert[4];

            gam1p := f(x0+h, V1p, Vp1p, get_v_delayed(x0+h-epsilon));

            gam2p := f(x0+2*h, V2p, Vp2p, get_v_delayed(x0+2*h-epsilon));

            g1p := g(x0+h, V1p, Vp1p, get_v_delayed(x0+h-epsilon),

                     (get_v_delayed(x0+h-epsilon+1e-8)-get_v_delayed(x0+h-epsilon))/1e-8);

            g2p := g(x0+2*h, V2p, Vp2p, get_v_delayed(x0+2*h-epsilon),

                     (get_v_delayed(x0+2*h-epsilon+1e-8)-get_v_delayed(x0+2*h-epsilon))/1e-8);

           

            F1p := h*vp0 - (V1p - v0 + h^2*(bet00*gam0 + bet10*gam1p + bet20*gam2p)

                   + h^3*(rho01*g0 + rho11*g1p + rho21*g2p));

            F2p := h*Vp1p - (V1p - v0 + h^2*(bet01*gam0 + bet11*gam1p + bet21*gam2p)

                   + h^3*(rho01*g0 + rho11*g1p + rho21*g2p));

            F3p := h*Vp2p - (V1p - v0 + h^2*(bet02*gam0 + bet12*gam1p + bet22*gam2p)

                   + h^3*(rho02*g0 + rho12*g1p + rho22*g2p));

            F4p := V2p - (2*V1p - v0 + h^2*(bet0*gam0 + bet1*gam1p + bet2*gam2p)

                   + h^3*(rho0*g0 + rho1*g1p + rho2*g2p));

           

            Fp := Vector([F1p, F2p, F3p, F4p]);

            J[1..4, j] := (Fp - F)/delta;

        end do;

       

        dV := LinearAlgebra:-LinearSolve(J, -F);

        V1 := V1 + dV[1]; V2 := V2 + dV[2]; Vp1 := Vp1 + dV[3]; Vp2 := Vp2 + dV[4];

    end do;

   

    return [V1, V2, Vp1, Vp2];

end proc:

 

# Main variable step-size loop

printf("Variable step-size integration for Example 1\n");

printf("tol = %e, h_init = %f\n", tol, h_init);

 

while x_curr < b - 1e-12 do

    # Try current step

    sol := solve_block(x_curr, v_curr, vp_curr, h_curr, omega);

    V1 := sol[1]; V2 := sol[2]; Vp1 := sol[3]; Vp2 := sol[4];

   

    # Compute with two half-steps

    sol_half1 := solve_block(x_curr, v_curr, vp_curr, h_curr/2, omega);

    V_mid := sol_half1[2]; Vp_mid := sol_half1[4];

    sol_half2 := solve_block(x_curr + h_curr/2, V_mid, Vp_mid, h_curr/2, omega);

    V2_half := sol_half2[2];

   

    # Error estimate

    err := abs(V2 - V2_half) / (2^6 - 1);

   

    if err < tol then

        # Accept step

        x_next := x_curr + 2*h_curr;

        X := [op(X), x_curr + h_curr, x_next];

        V := [op(V), V1, V2];

        Vp := [op(Vp), Vp1, Vp2];

        x_curr := x_next;

        v_curr := V2;

        vp_curr := Vp2;

        printf("x = %7.4f, h = %8.5f, err = %12.5e\n", x_curr, h_curr, err);

       

        # Adjust step size

        if err < tol/2 then

            h_curr := min(2*h_curr, h_max);

        end if;

    else

        # Reject step, reduce h

        h_curr := max(h_curr/2, h_min);

        printf("  Rejecting, new h = %8.5f\n", h_curr);

    end if;

end do:

 

# Exact solution for comparison

exact := x -> 3*sin(x) - 5*cos(x);

errors := [seq(abs(V[i] - exact(X[i])), i=1..nops(X))];

 

# Visualization

p1 := pointplot([seq([X[i], errors[i]], i=1..nops(X))], color=red, symbol=circle,

                title="Example 1: Variable Step-Size TFSB - Absolute Errors",

                labels=["x", "Error"], labeldirections=[horizontal,vertical]);

p2 := plot([[x_curr, h_curr]], x=a..b, style=point, color=blue,

            title="Step-size evolution", labels=["x", "h"]);

display(p1);

display(p2);

 

printf("\nFinal results for Example 1:\n");

printf("Number of steps: %d\n", nops(X)-1);

printf("Maximum error: %e\n", max(errors));

printf("Final step-size: %f\n", h_curr);

Since Maple version 2026, I have noticed that the orientation option with the contourplot3d command has no effect.

It is easy to control by using the example provided in the online help and adding, for example, orientation=[20,10,10].

Thank you for your help.

Best regards.

On my system (Windows 11) I can only remove an entry from the favorites when a document is open. With all documents closed removing does not work.

Can someone confirm?

Color option in Context Panel does not work

Steps to Reproduce:
1. Create or open a histogram.
2. Try to change its color using the Color option in the Context Panel.
3. Observe that the color does not change.
4. Click outside the histogram.
5. Click back on the histogram.
6. Try the Color option again — now it works.

Expected Behavior:  
The Color option should work immediately when applied to the histogram, without requiring extra clicks.

Actual Behavior:  
The Color option only works after clicking outside the histogram and then reselecting it.

Just for my interest.
Why is the following not working

one := ``(1);
                           one := (1)

 lprint(`%`);
Error, Got internal error in Typesetting:-Parse:-Preprocess : "invalid subscript selector"
Typesetting:-mambiguous(Typesetting:-mambiguous( lprint(%), 

  Typesetting:-merror("Got internal error in Typesetting:-Parse:\

  -Preprocess : "invalid subscript selector"")))

but this works

lprint(one);
``(1)

I do not understand why a once entered prompt cannot be copied from the AI pannel.

I also do not understand that a text passage from an answer (e.g. a proposed help topic) cannot be selected and copied to the clipboard with crtl-c.

It does not make sense to me to unnecessarily restrict a new program feature. This somehow spoils the show.

Is that a restriction to Windows or my local setting or something imposed on Maplesoft from a third party supplier?

Has anybody seen similar AI implementations in other products?

Hello Ladies and Gents :)

I asked Maple AI to write a program to calculate the geodesics on a sphere.

see the attached ws.

When running this program I receive an error which I don't understand.

of course i know how to do it by paper and pen, it was just a test.

if someone could explain me this error i would greatly appreciate .

thank you and kind regards,

Jean-Michel

Moderator edit - here is the file:

geodesics.mw

Here is another case where an incorrect solution of a non-algebraic expression is returned.

How to make Maple exclude solutions obtained with allsolutions where the expression is not defined?

IMO, allsolutions facilities of Maple should do this automatically. Anything from a mathematical point of view that speaks against an automatic exclusion?

(Are there other commands that provide a correct solution?)

By the way:

solve(sin(x)/x)

does not return a single solution as stated on the solve help page ("In general for transcendental equations, the solve command returns only one solution,...") because solve correctly discards x=0 but does not consider returning x=pi.

>  >  (1) The above expression is not defined at >  (2) >  Error, (in value/Eval) numeric exception: division by zero   Roots of the expression with RootOf >  (3) >  (4) >  Originally _Z1, renamed _Z1~:   is assumed to be: integer     >  (5) >  (6) >  (7) >  (8) The above should be false because  is not defined for However, this is what automatic simplifcation does behind the scenes with the output (7)  >  (9) For this particular case Maple should have returned a special name expressing nonzero integers, like  . or or are common symbols for that So is incorrect in the solution and automatic simplification does something wrong. Only substituing everything at once leads to a correct response: >  Error, numeric exception: division by zero   Solve >  >  (10) >  (11) >  (12) >  (13) >  (14) >  Error, numeric exception: division by zero   Again is not a valid solution for the first solution butworks for the second solution >  (15) >  (16) >  (17) Substituing everything at once works >  (18) >  (19) However, the second solution does not cover all solutions (it misses even multiples of π) >  (20) >  (21)

Related topic: I also think that new users and students should profit from a directly understandable output like this

Maybe adding new convert form is an option.

Download roots_of_sinc.mw

I do have a stundent license permanently and I've installed the OpenMaple java examples: JCEngineCallBacks, JCHelpCallBack, jcmaple in Eclipse.

I've also installed the MapleTest.java which runs fine.

When I start the jcmaple.java I get the error: license communication problem 1  Error starting OpenMaple session

So I don't know why.

Meanwhile I was able to solve the issue

...and starting values. I want to reproduce the maximum in Maple using the solution structure in the attached file, similar to the old Mathcad method. Instead of using the Euler-Lagrange equation, I want to separately enter the objective function, then the boundary conditions, and finally the starting values ​​for an iteration (as was done in the old Mathcad solution block). I entered the latter into the "maximize" command using the Maple help text. But nothing is being calculated. What am I doing wrong? As I said, I only want a numerically generated approximate solution.test.mw

>  >  (1) >  (2) >  (3) >  (4) >    >  >  (5) >  (6) >  (7) >  (8) >  (9)   >  (10)

Download test.mw

The command completion facility is real time saver. On the other hand, entering the colon dash is not fluent since it requires searching the keys on the key board.

I was wondering if there is a keyboard shortcut that enters both characters at once.

Is there an easier or shorter way to do the following? Any handy package command?

(1) (2) (3) (4)

Related question: Is that substitution mathematically correct or does the _a belong to the outermost RootOf?

Download parameters_in_nested_RootOf.mw