What you experience are numerical integrations errors.

method=rosenbrock gives a slightly deceasing energy.

Forward Euler is in this respect particulary bad (and maybe instructive for your students). I came accross a similar behaviour when exloring solver settings in MapleSim

In this context I ask myself: What has to be invested in numerical fidelity when it comes to predictions of the future of the solar system?

Below is way using names p1 and m1 for 1 and -1 to prevent automatic simplification. These names are assigned to other names that typeset nicely. This assignement is not mandatory but facilitates input (macro could have been used alternatively).

In addition auxilliary equation expressions are used to counteract automatic simplification (ax1) and to simplify names with brackets (ax2). The equation expressions are generally valid and do not create overhead. All other manipulations are performend on equations (i.e. respecting equality).

(Initially I tried to work with inert operators but this did not work well because I could not expand such expressions.)

Download m1_times_m1.mw

What about:

e:=_C1^2+_C1+_C2^3+a+b;
indets(e, And(symbol, suffixed(_C, nonnegint))^anything);#extract powers of
indets(e-add(%), And(symbol, suffixed(_C, nonnegint))); #remove powers of and scan for remainders
%% union %

{_C1, _C1^2, _C2^3}

Edit:
This one works also for products

e:=_C1^2+_C1+_C2^3+a+b:
pw_ind:=indets(e, And(symbol, suffixed(_C, nonnegint))^anything):
lin_ind:=indets(subs(seq(pw_ind[i]=0,i=1..nops(pw_ind)),e), And(symbol, suffixed(_C, nonnegint))):
pw_ind union lin_ind
                       /        2     3\ 
                      { _C1, _C1 , _C2  }
                       \               / 

Maple assumes all names complex. Just look at the two first terms of the sum when we force Re on it

[op(1..2,lhs(P))];
map(Re,%);
add(%);

Without assumptions Maple cannot tell which term is real and which is imaginary.

PS.: In one of the terms there might be a space missings

In case you are interested in a way to get all solutions in one go (Roubens answer contains them in the text passage: i.e.: m=0 is allways a solution and the two other solutions are only opposite in sign)

The variant below creates a list of time stamps, solves your equation eqm for m, substitutes both into eqG and sorts the results.

eqm := m = tanh(6*m/T);
eqG:=G = -T*ln(2.0*cosh(6*m/T)) + 3*m^2;# 
T_list:=[seq(t,t=0.1..4,0.1)]:
m_list:=[seq((op~)([fsolve(subs(T=t,eqm),{m=-2..2},maxsols =3)]),t in T_list)]:
G_list:=[seq([seq(evalf(subs(m_list[j][i],T=T_list[j],eqG)),i=1..nops(m_list[j]))],j=1..nops(T_list))]:
for i to nops(T_list) do T=T_list[i],zip((x,y)->[x,y],m_list[i],G_list[i])[]; end do;

I have derived the answer from an attempt of using solve and getting all the roots from a RootOf expression of the symbolic solution. Neither allvalues nor evalf(RootOf) worked as I had hoped-for. In the attached document I have solved the argument of the RootOf expression with fsolve. This turned out the be numerically more demanding (Digits had to be increasesd to find all roots) than directly solving the problem numerically (as above). Such a hybrid approach might be usefull in other situations but not here.

tanh_cx_-_x.mw

Using op, you can disasemble expressions into their operands. I have highlighted this in yellow below.

The result represents a sum over _alpha where the summation index _alpha equals a RootOf expression. The RootOf expression is a placeholder for all roots of an expression. For some expressions RootOf can be evaluated.See ?RootOf for more. In your case you get 4 roots to choose from. 

allvalues(RootOf(1 + C1*C2*L1*L2*_Z^4 + (C1*C2*L1*R2 + C1*C2*L2*R1)*_Z^3 + (C1*C2*M^2*omega^2 + C1*C2*R1*R2 - C1*L1 - C2*L2)*_Z^2 + (-C1*R1 - C2*R2)*_Z));

This answer comes delayed (MaplePrimes service was not available for me).

If op is an option for you:

op([1,1],Zin)*op(2,Zin)+op([1,2],Zin)*op(2,Zin)

or a bit more generic for longer sums (if map is acceptable as well)

map(op(0,Zin), op(1,Zin), op(2,Zin))

For code readability and avoidance of errors, I am not a big fan of using op. In this case it might be an acceptable compromise. Whenever I have to use op, I check the simplified result for equivalence.

acers variants shine a light on how a custom-made "simplification" routine could be implemented not using "low level" command as op.

I agree that a command would be desirable that expands a quotient (having a sum in its numerator) without applying any further algebraic transformation to the expression (as denom and numer do to Zin).

I would also welcome a command that can collect more complex expressions (like the denominator in your example) from an expression than it is possible today.

The dot needs to be removed. Then you will see that d is a table

with two entries that should be the same but are not. The assignement above is problematic. If you assign d[1] and d[2] to something else then the mysterious t in the output of the evalc call is gone.

The error message indicates that pdsolve cannot evaluate the initial condition with the integral. Here is a way to make your code  work

restart;
with(PDEtools);

# the heat equation
pde := diff(v(x, t), t) = diff(v(x, t), x, x);

f := x -> exp(-x^2);

a := -1.68858; b := 1.68858; delx := 1/100; t_final := 5;

pds := pdsolve(pde, {v(a, t) = 0, v(b, t) = 0, v(x, 0) = (1 - convert(series(int(f(f(s)), s = 0 .. x),x),polynom))}, numeric, range = a .. b, time = t, spacestep = delx);

pds:-animate(frames = 120, t = t_final);

This techinque is an option in cases where analytic solutions are not availabe or initial conditions cannot be evaluated. The order of the series approximation determines (of course) the fidelity of the solution - worth comparing to acers solution which is more precise.

You have declared some global variables as local. To evaluate expressions with gobal variables with values of local variables do

sysrr := eval(totalsecular3, [:-Q__1 = Q__1, :-Q__2 = Q__2, :-Omega__1 = Omega__1, :-Omega__2 = Omega__2, sigma = varr1, r = 0.047])

Use the return command or the print command to "output" Q1 and Q2.

You still need to fix an error with the Jacobian. I leave this up to you since I do not understand what you mean by "stability of judgement".

First example: Maple seems to apply a limit to negative infinity to a solution with a finite boundary which leads to y(x)=0.

With a finite boundary the solution seems to match the numerical solution from Paras31.
Also, the symbolic solution for a finite boundary passes odetest.

IMO, Maple should not have returned this solution for the boundary values because the slope is undefined at infinity for the general solution of the ode.

diff(dsolve(ode),x):
eval(%,x=-infinity)
              / d                       \            
          eval|--- y(x), {x = -infinity}| = undefined
              \ dx                      /            

wrong_sol_dsolve_sept_10_2024_with_limiit.mw

Frequencies are correct but the plots are against expectation at the first glance.

The oscillation period for the simulations at the length of 0.001 and 1.001 are the same when you zoom in your plots. That is what can be expected for the pendulum and better be seen on the plot of Rouben Rostamian. Compare for example in your worksheet

odeplot(sol1, [t, theta(t)], 3998 .. 4000);
odeplot(sol2, [t, theta(t)], 0. .. 2);

What is not intuitive at the first glance is the change in amplitude:

Lengthening the pendulum (the red curve) increases the angular excursion.

Update: I deleted the explanation for the increase of amplitude because the ode does not correctly describe a pendulum of variable length. A simulation with MapleSim of a pendulum (pointmass attached to a string of variable length) does not show the increase in amplitude. (However, the frequencies are consistent with odeplot results.)

The ode of the original post was derived for a perfect circular motion. It cannot describe the dynamics of a lengthening pendulum. Introducing a variable length after assuming a constant length does not make sense. For a slow change of the length osciallation frequencies are well described but the magnitude of the angular excursion must be ignored.

You could use this alternatively

Optimization:-Minimize(G(x, y), x = 0 .. 1, y = 0 .. 1);
       [-0.710208044926744719, 

         [x = 0.894924318449926, y = 0.971579686253484]]

For the functions you have to load a Student package. Paste this in an input line

with(Student[VectorCalculus]);
Hessian(sqrt(81 + (u + 5)^2) + sqrt(49 + (v - 3)^2) + sqrt(4 + (v - u + 1)^2), [u, v]);
Gradient(sqrt(81 + (u + 5)^2) + sqrt(49 + (v - 3)^2) + sqrt(4 + (v - u + 1)^2), [u, v]);

For the long brackets and > push on or Insert -> Execution group

The semicolon terminates a statement in an input line. It's called the statement separator. Enter help(";") for more.

For a 3d plot enter ?plot3d

Before doing that you might want to follow the startup page. 

If you cannot find it, navigate to View->Home

Restart resets Maple. It basically "forgets" all input and computation results.

According to equation (13) of ?odeadvisor,Chini the following ode is Chini

ode := diff(y(x),x)=a*y(x)^n+b*x^(n/(1-n));
DEtools:-odeadvisor(ode)
                                            /  n  \
                                            |-----|
                     d               n      \1 - n/
             ode := --- y(x) = a y(x)  + b x       
                     dx                            

               [[_homogeneous, class G], _Chini]

Meaning for me in

Chini_ode := diff(y(x), x) = f(x)*y(x)^n - g(x)*y(x) + h(x)
                   d                  n                   
     Chini_ode := --- y(x) = f(x) y(x)  - g(x) y(x) + h(x)
                   dx                                     

 that f, g , h can be constant (functions)
 

f(x)=a,g(x)=0,h(x)=b*x^(n/(1 - n));
subs(%,Chini_ode);
                                           /  n  \
                                           |-----|
                                           \1 - n/
             f(x) = a, g(x) = 0, h(x) = b x       

                                        /  n  \
                                        |-----|
                 d               n      \1 - n/
                --- y(x) = a y(x)  + b x       
                 dx                            

My understanding from the explanations from @ecterrab is that odeadvisor does not necessarily returns all matching types. In your case Chini is maybe not considered an appropriate method and therefore not suggested. This would fit to your observation of a division by zero error.

What is still unclear for me is the otherwise in: "odeadvisor might return more than one type; otherwise, the first matching of a pattern interrupts the process and a classification is returned". Why is in this case the pattern matching process interrupted in others (see above) not.

`dsolve/methods`[1]
  [quadrature, linear, Bernoulli, separable, inverse_linear, 

    homogeneous, Chini, lin_sym, exact, Abel, pot_sym]

If called without options, odeadvisor is not a classification command that retruns all types. It is supposed to be a "classification and suggestion" command. Meaning it retruns all classifications of suggested methods to solve. If this is the correct interpretation, the help page needs a bit more clarification. I am saying that because I read it serveral times and interpreted it differently in the beginning.