dharr

AIC...

Answered: dharr 8812

December 17 2025

1 3

I didn't know about symbolic regression, but it seems from your responses to @sand15 that you might also be interested in automatic selection from a class of models. He has provided a nice summary of the issues. Since he mentioned AIC and provided nice code for rational function fitting, I thought I would illustrate the principle of AIC for those functions.

Basically the idea is that more fitting parameters generally improve the fit but after a while the extra parameters are not meaningful (overfitting). AIC quantifies this tradeoff into one parameter, AIC. I recklessly ignored the full rank warnings; this is just to give you some idea of the method.

One can also use an F-ratio test to see whether a model with another parameter is statistically significant. Here you need to choose a confidence level, and it is not always clear what to use. I used to use F-ratio (at the 1% level), but more recently switched to AIC, though I have the impression that it still overfits more that I would like; that may just prove I'm too cautious. AIC can select between different classes of models, but you have to generate them in some way, which I guess is what symbolic regression does.

At the end of the day your data looks suspiciously noise-free, and the key issue may be more of finding what type of function would look like that (with less that the 16 parameters AIC decided on) than the regression itself. That still seems like a human endeavour.

@sand15's FitRationalFraction procedure and the dataset are in the startup code

>

restart;

>	kernelopts(version);

>	X := <map2(op, 1, dataset)>: Y := <map2(op, 2, dataset)>: n := nops(dataset);

Fit all rational fraction models with numerator/denominator degrees between 5 and 10

>	AIC:=Array(5..10,5..10):

>	for ij in indices(AIC) do sij := ij[]; F := FitRationalFraction(ij, X, Y, AllResults=true): S := F[1]["residualsumofsquares"]; p := numelems(F[1]["parametervalues"]); fit[sij] := F[2]; AIC[sij]:=evalf(nln(2Pi)-nln(n)+n+nln(S)+2*(p+1)); end do:

Warning, model is not of full rank

In all rows and columns, the AIC gets better (more negative) and then worse (as overfitting occurs).

>

AIC;

`Array(5..10, 5..10, {})`

We want the fit with the most negative AIC

>	min(AIC); min[index](AIC); bestfit:=fit[%];

>	plots:-display( plot(rhs(bestfit), x=(min..max)(X), color="SteelBlue") , Statistics:-ScatterPlot(X, Y, symbol=circle, symbolsize=10, color=black) , view=[(min..max)(X), (min..max)(Y)] , size=[1000, 400] )

>

Download RationalFractionFitAIC.mw

A way...

Answered: dharr 8812

December 12 2025

2 0

I didn't try to figure out which solutions you might want; most of them have many a[i] and b[i] zero, so you may want to solve for different variables.

>

(1)

>

ode := (diff(diff(diff(diff(U(xi), xi), xi), xi), xi))*`ε`-(-6*w^2*`ε`+alpha)*(diff(diff(U(xi), xi), xi))-U(xi)*(3*w^4*`ε`-alpha*w^2-4*mu*U(xi)^2+k)

(2)

>

F := sum(a[i]*G(xi)^i, i = 0 .. 2)+sum(b[i]*G(xi)^(-i), i = 1 .. 2)

a[0]+a[1]*G(xi)+a[2]*G(xi)^2+b[1]/G(xi)+b[2]/G(xi)^2

(3)

>

(4)

>

(5)

Look, ma, no square roots!

>

U(xi) = (a[2]*G(xi)^4+a[1]*G(xi)^3+a[0]*G(xi)^2+b[1]*G(xi)+b[2])/G(xi)^2

(6)

>

(7)

Still have derivatives of G, so

>

(8)

>

(9)

>

${4*mu*b[2]^3, 12*mu*b[1]*b[2]^2, 120*epsilon*h^4*a[2]+4*mu*a[2]^3, 12*mu*a[0]*b[2]^2+12*mu*b[1]^2*b[2], 24*epsilon*h^4*a[1]+12*mu*a[1]*a[2]^2, 24*mu*a[0]*b[1]*b[2]+12*mu*a[1]*b[2]^2+4*mu*b[1]^3, -120*epsilon*h^4*a[2]-36*epsilon*h^2*w^2*a[2]+6*alpha*h^2*a[2]+12*mu*a[0]*a[2]^2+12*mu*a[1]^2*a[2], -20*epsilon*h^4*a[1]-12*epsilon*h^2*w^2*a[1]+2*alpha*h^2*a[1]+24*mu*a[0]*a[1]*a[2]+4*mu*a[1]^3+12*mu*a[2]^2*b[1], 16*epsilon*h^4*b[2]+24*epsilon*h^2*w^2*b[2]-3*epsilon*w^4*b[2]-4*alpha*h^2*b[2]+alpha*w^2*b[2]+12*mu*a[0]^2*b[2]+12*mu*a[0]*b[1]^2+24*mu*a[1]*b[1]*b[2]+12*mu*a[2]*b[2]^2-k*b[2], epsilon*h^4*b[1]+6*epsilon*h^2*w^2*b[1]-3*epsilon*w^4*b[1]-alpha*h^2*b[1]+alpha*w^2*b[1]+12*mu*a[0]^2*b[1]+24*mu*a[0]*a[1]*b[2]+12*mu*a[1]*b[1]^2+24*mu*a[2]*b[1]*b[2]-k*b[1], 16*epsilon*h^4*a[2]+24*epsilon*h^2*w^2*a[2]-3*epsilon*w^4*a[2]-4*alpha*h^2*a[2]+alpha*w^2*a[2]+12*mu*a[0]^2*a[2]+12*mu*a[0]*a[1]^2+24*mu*a[1]*a[2]*b[1]+12*mu*a[2]^2*b[2]-k*a[2], epsilon*h^4*a[1]+6*epsilon*h^2*w^2*a[1]-3*epsilon*w^4*a[1]-alpha*h^2*a[1]+alpha*w^2*a[1]+12*mu*a[0]^2*a[1]+24*mu*a[0]*a[2]*b[1]+12*mu*a[1]^2*b[1]+24*mu*a[1]*a[2]*b[2]-k*a[1], -8*epsilon*h^4*b[2]-12*epsilon*h^2*w^2*b[2]-3*epsilon*w^4*a[0]+2*alpha*h^2*b[2]+alpha*w^2*a[0]+4*mu*a[0]^3+24*mu*a[0]*a[1]*b[1]+24*mu*a[0]*a[2]*b[2]+12*mu*a[1]^2*b[2]+12*mu*a[2]*b[1]^2-k*a[0]}$

(10)

>

(11)

>

(12)

>

(13)

>

(diff(diff(diff(diff(U(xi), xi), xi), xi), xi))*epsilon-(-6*w^2*epsilon+(epsilon*h^4+6*epsilon*h^2*w^2-3*epsilon*w^4-k)/((h-w)*(h+w)))*(diff(diff(U(xi), xi), xi))-U(xi)*(3*w^4*epsilon-(epsilon*h^4+6*epsilon*h^2*w^2-3*epsilon*w^4-k)*w^2/((h-w)*(h+w))+k)

(14)

>

U(xi) = b[1]/sech(h*xi)

(15)

>

(16)

Download test-F-p2.mw

simpler...

Answered: dharr 8812

December 09 2025

1 3

@sand15 I think you meant rhs(sol) in unapply. Then the second step can be simpler:

Sol := unapply(rhs(sol), x);
eval(IC, y = Sol);

Collecting with symbolic exponents...

Answered: dharr 8812

December 07 2025

4 0

Edit: handling of terms like ln(r) improved.

https://www.mapleprimes.com/questions/242062-Collect-Using-Pattern

>

restart;

We want to collect A with respect to r, including symbolic exponents, in order to get B. I've added a constant term to the OP's example, to further test the method.
It is not clear from the question if the order matters, or what it would be for more complicated exponents.

For simplicity, assume that the order of exponents is just the default sort order, which would work here.

>	A:=r^(2a)+r^2(1+a+r^(2a)) + r + ar + a; B:=a+(1+a)r+(1+a)r^2+r^(2a)+r^(2+2a);

Notice that Maple's automatic simplification may frustrate any attempt to get a desired order, e.g., the following is already collected, but does not retain the input order, and is perhaps not in the expected order

>	r^3 + a*r^2 + r;

The reason collect does not work is that it is for polynomials, which have integer exponents but not symbolic ones.
The following finds the exponent wrt r for any exponent not containing r

>	ex := (y,r) -> r*diff(y,r)/y;

proc (y, r) options operator, arrow; r*(diff(y, r))/y end proc

Examples. The presence of r in the second case indicates we should probably just leave this alone (this is done in the procedure below).

>	ex(7r^(2tan(a)), r); ex(a*ln(r), r);

Expand, find the terms and their corresponding exponents and coefficients. (assumes it is type `+`)

>	terms := [op(expand(A))]; exps := expand(ex~(terms, r)); cfs := zip((t,e)->expand(t/r^e), terms, exps);

Sort them

>	exps2, p := sort(exps, output = [sorted, permutation]); cfs2 := cfs[p];

Assemble the results

>	s:=0: cf:=cfs2[1]: expt:=exps2[1]: for i from 2 to nops(exps2) do if exps2[i] = expt then cf += cfs2[i] else s += cfr^expt; expt := exps2[i]; cf := cfs2[i] end if; end do; s += cfr^expt;

As a procedure.

>

Collect := proc(expr, x::name)
local zz, ex, term, cfs, exps, p, s, cf, i;
if expr::function then return expr end if;
for i,term in expand(expr)+zz do # zz forces type `+`
  ex := expand(x*diff(term,x)/term);
  if has(ex, x) then
     cfs[i] := term;
     exps[i] := 0;
  else
     exps[i] := ex;
     cfs[i] := expand(term/x^ex)
  end if;
end do;
exps := convert(exps, list);
cfs := convert(cfs, list);
exps, p := sort(exps, 'output' = ['sorted', 'permutation']);
cfs := cfs[p];
s := 0;
cf := cfs[1];
ex := exps[1]:
for i from 2 to nops(exps) do
  if exps[i] = ex then
    cf += cfs[i]
  else
    s += cf*x^ex;
    ex := exps[i];
    cf := cfs[i]
  end if;
end do;
s += cf*x^ex - zz;
end proc:

Test

>	Collect(A + ln(r), r);

>	Collect(r^3 + a*r^2 + r, r);

>	Collect(ln(x)+3*x^2, x);

x^f(x) is left alone, i.e., is treated as a constant term.

>	Collect(3x^2 + 2x^x, x);

>

Download CollectSymbolicExponents.mw

simplify with side relations...

Answered: dharr 8812

December 04 2025

1 3

sqrt(1-cos(x)^2);
simplify(%,{cos(x)^2=1-sin(x)^2});

expand...

Answered: dharr 8812

December 01 2025

1 0

For some reason, the original simplify gives a numerator/denominator, with sqrt(1-t^2) in both. If expand is used first then we don't have numerator/denominator, Then simplify can deal with the diffs. As well the denominator sqrt(1-t^2) might cancel or be simpler. So simplify(expand(B)) works here. The collect serves the same purpose as expand - giving priority to the diffs means we get a polynomial in diffs and not numerator/denominator, which then simplifies. The collect can be just wrt diff, so simplify(collect(B,diff))works.

Error in paper...

Answered: dharr 8812

November 23 2025

0 1

As far as I can see, there is some error in the paper. Maybe one of the other case(s) works.

Download pde.mw

Numbers...

Answered: dharr 8812

November 19 2025

3 1

For the numbers.

Download Numbers.mw

For the calling of previous commands by referring to them by the equation labels on the right-hand side, use ctrl-L to insert the number in the appropriate place.

open panel/palettes...

Answered: dharr 8812

November 18 2025

3 1

If you seach for "slow" on this site, there are mentions that having the context panel - the panel on the right - closed can increase performance; you might also try closing the palettes on the left when you are not using them.

Elliptic curve solutions...

Answered: dharr 8812

November 16 2025

4 0

[Edited]

parametrization is good for finding the first rational point (not necessarily integer) on genus 0, but doesn't work for genus 1.

I assumed at first you wanted integer solutions. My usual attacks here were not that successful. For the first equation 1 found 1 more solution, but rather trivial, and then the addition process gives the infinity point. For the second case, it seems to be working nicely in producing rational solutions.

Elliptic.mw (first equation)

Second equation:

https://www.mapleprimes.com/questions/242001-Elliptic-Curves

>

Find integer solutions of the following equation of interest; we know one solution is (3,7);

>

Check genus is 1

>

Define addition of points in the group algebra - see Wikipedia

A is the coefficient in y^2=x^3+A*x+B

>

algadd := proc(P::[rational,{rational,identical(infinity)}], Q::[rational,{rational,identical(infinity)}], A::rational)
# infinity point denoted [0, infinity]
  local s,xP,yP,xQ,yQ,xR,yR;
  xP,yP := P[];
  xQ,yQ := Q[];
  if yQ = infinity and yQ = infinity then return [0, infinity] end if;
  if yP = infinity then return [xQ, yQ] end if; # infinity point is identity
  if yQ = infinity then return [xP, yP] end if;
  if xP <> xQ then
    s := (yP - yQ)/(xP - xQ);
    xR := s^2 - xP - xQ;
    yR := yP - s*(xP - xR);
    [xR, -yR];
  elif yP = -yQ then
    [0, infinity] # infinity point
  else # P = Q,sum is tangent
    s := (3*xP^2 + A)/(2*yP);
    xR := s^2 - 2*xP;
    yR := yP - s*(xP - xR);
    [xR, -yR]
  end if;
end proc:

The j_invariant has the same value despite the various transformations (birational equivalence)

>

Maple finds a Weierstrass form

>

wform := Weierstrassform(f, x, y, X__1, Y__1); w1 := wform[1]; j_invariant(w1, X__1, Y__1)

The addition routine wants it in the standard form . For integer solutions, take A and B integers, which we can achieve by changing the sign of X, scaling the whole thing by 9*81 and redefine X and Y by factor of 9 and 27. (For rational solutions, rational A and B are OK, so we only need to change the sign of X here.)

>	$trans := {X__1 = -(1/9)X, Y__1 = (1/27)Y}; w := eval((981)w1, trans); j_invariant(w, X, Y)$

${X__1 = -(1/9)*X, Y__1 = (1/27)*Y}$

Conversion between f and w

>	$itrans := solve(trans, {X, Y}); xy_to_XY := eval(itrans, {X__1 = wform[2], Y__1 = wform[3]})$

${X = -9*X__1, Y = 27*Y__1}$

Initial solution in X,Y

>

Conversion between w and f

>	$XY_to_xy := eval({x = wform[4], y = wform[5]}, trans); eval(XY_to_xy, sol1XY)$

To start the algorithm, we need A and the first solution

>

isolate(w, Y^2); A := coeff(rhs(%), X, 1); P[1] := eval([X, Y], sol1XY)

Add P[1] to get another solution - it's rational but not integer

>

P[2] := algadd(P[1], P[1], A); eval(XY_to_xy, `~`[`=`]([X, Y], %)); eval(f, %)

>

P[3] := algadd(P[2], P[1], A); eval(XY_to_xy, `~`[`=`]([X, Y], %)); eval(f, %)

>

P[4] := algadd(P[3], P[1], A); eval(XY_to_xy, `~`[`=`]([X, Y], %)); eval(f, %)

To find the first solution, we want either , or a Y value for which divides the disciminant, for which the corresponding X is an integer (Nagell-Lutz) .

>

Cubic := rhs(isolate(w, Y^2)); d := discrim(Cubic, X); ifactor(d)

So test 0 and all the divisors of 2*3^6

>	$divs := `union`({0}, NumberTheory:-Divisors(mul(map(proc (q) options operator, arrow; ifelse((q[2])::even, q[1]^((1/2)*q[2]), 1) end proc, ifactors(d)[2]))))$

>

j := 0; for i in divs do facs := factors(eval(w, Y = i))[2]; if nops(facs) > 1 then j := j+1; sol[j] := {X = solve(select(proc (z) options operator, arrow; type(z[1], linear) end proc, facs)[1][1], X), Y = i}; print(sol[j]) end if end do

No luck here.

>

Download Elliptic2.mw

timelimit...

Answered: dharr 8812

November 15 2025

2 10

A few of points:

1. Sometimes a calulation can be fast but displaying the complicated expression can be slow, so it is better then to store the answer without displaying it, e.g. use ans := pdetest(eqv2, pde): then if it completes you can attempt to display it (perhaps check its length first).

2. You can use ans := timelimit(20, pdetest(eqv2, pde)): to limit the attempt to 20 s.

3. If you suspect a complicated expression might simplify to zero you can set random numerical values for parameters and see if the numerical value is nearly zero. In your case you have functions beta(t) etc, so you might give simple functions to these, e.g.,

params := {beta=(t->t^2), k[1] = 0.35, ...};
ans := eval(eval(pde, params), eval(eqv2, params)):

Converting to a first order system...

Answered: dharr 8812

November 14 2025

0 3

I didn't undertand what you said. It looks as though you were trying to manually make a first order system. I fixed this and get the same result as convertsys.

systems.mw

bad parameter values...

Answered: dharr 8812

November 13 2025

2 5

Download find-parameter.mw

Bug...

Answered: dharr 8812

November 12 2025

1 0

Someone else found a similar problem here and I reported this bug to Maple. The same workaround might work here - put simplify(...) around the expression with the square root.