tomleslie

Tidied up a few things...

Answered: tomleslie 13876

February 21 2019

0 0

The 'if' test in Fract contains a redundant 'is()' - if condition are treated as Boolean
Fixed the equation in Fract1
Handled the issue in Fract1 where isolve() can return an infinite number of solutions
Added a 'No Solution' clause where appropriate
In the attached, both procedures give the same answer

>

  Fract:= proc( P::posint, Q::posint)
                local p, q;
                for p to P-1 do
                    for q to Q-1 do
                        if   (P-p)*q-p*(Q-q) = 1
                        then return p/q, P/Q, (P-p)/(Q-q)
                        end if
                    end do:
                end do:
                return `No Solution`
          end proc:
  Fract1:= proc( P::posint, Q::posint)
                 local p, q, sol;
                 sol:= eval( isolve((P-p)*q-p*(Q-q) = 1),
                             _Z1=0
                           );
                 if   sol <> NULL
                 then return eval( [p/q, P/Q, (P-p)/(Q-q)],
                                   sol
                                 )[];
                 else return `No Solution`
                 fi;
           end proc:
  Fract(7, 81);    Fract1(7,81);
  Fract(39, 97);   Fract1(39,97);
  Fract(101, 143); Fract1(101,143);
  Fract(11, 80);   Fract1(11,80);
  Fract(15, 37);   Fract1(15,37);
  Fract(22, 39);   Fract1(22,39);
  Fract(25, 37);   Fract1(25,37);
  Fract(21, 91);   Fract1(21,91);
  Fract(13, 19);   Fract1(13,19);

(1)

>

Download dioph.mw

Makes no difference...

Answered: tomleslie 13876

February 19 2019

1 0

If you use

showstat(DocumentTools:-GetProperty);

to examine the code for the GetProperty() command, the first line is

DocumentTools:-GetProperty := proc(id::{string, symbol}, prop::{string, symbol, typeindex(symbol)})

so the type of the second argument has to be either a string a symbol, or an indexed symbol it doessn't matter which. So in your case either value or "value" is fine

Well...

Answered: tomleslie 13876

February 19 2019

1 0

It may depend on whther you are going to "generalize" this code in any way. For the case you provide, the following provides the same answer (and is shorter)

Download someCombs.mw

For the example you give...

Answered: tomleslie 13876

February 19 2019

0 8

it depends whether you want an exact (ie symbolic) solution, or whether a numeric approximation is acceptable?

The attached finds both

>

  fapp:= x^(2)+423*x^(3)/3+6789*x^(4)/4;
  fexact:=sin(x)/x;
#
# maximize() will not find an a symbolic
# solution if the abs() function is used,
# but one can find the maximum square error
#
  maximize((fapp-fexact)^2, x=0..2, location)[2][];
#
# Need to take the square root to get back to
# the maximum absolute error
#
  sqrt(%[2]);
#
# Generate a floating point approximation for
# the above
#
  evalf(%);

$[{x = 2}, (1/64)*(-226304+4*sin(2))^2]$

(1)

>	# # Use the optimization package to find the maximum absolute error # numerically # Optimization:-NLPSolve( abs(fapp-fexact), x=0..2, maximize=true);

(2)

>

Download maxErr.mw

Don't cut/paste from worksheet...

Answered: tomleslie 13876

February 19 2019

0 2

Doing this just produces a garbled mixture of Maple input and Maple output, which is almost impossible to interpret.

Instead use the big green up-arrow in the Mapleprimes toolbar to upload your worksheet.

If I try tidying up the cr*p you did supply then I get the code shown in the attached. However your original post also contains references to a procedure given only as

rang := proc(M::list, a) ... end;;

where the body of the procedure is not supplied, and hence cannot be debugged by anyone here

>

restart;
A002487:= proc(m)
               local a, b, n;
               option remember;
               a:= 1; b:= 0; n:= m;
               while 0 < n do
                     if   type(n, odd)
                     then b:= a+b
                     else a:= a+b
                     end if;
                     n:= floor((1/2)*n)
               end do;
               b
          end proc:

listeinverse:= proc(L::list)
                    local i;
                    [seq(op(nops(L)-i, L), i = 0 .. nops(L)-1)]
               end proc:

Brocot:= proc(n)
              local c, i, L, M, r;
              L:= NULL;
              r:= 2^n;
              L:= [seq(A002487(i), i = 0 .. r)];
              M:= listeinverse(L);
              c[0]:= 0, 1/cat(0);
              for i to r do
                  c[i]:= L[i]/M[i]
              end do;
              c[r+1]:= 1/cat(0);
              return [seq(c[i], i = 1 .. r+1)], r+1
        end proc:

F:= proc(N)
         local a, b, L; L:= NULL;
         L:= sort([op({seq(seq(a/b, a = 0 .. b), b = 1 .. N)})]);
         return L, nops(L)
    end proc:

for i from 0 to 4 do
    B || i:= Brocot(i)
end do;

F(1); F(2); F(3); F(4);

(1)

>

Download newProc.mw

When I execute...

Answered: tomleslie 13876

February 18 2019

0 5

the worksheet Fuzzy_clusters_subtour.mw it generates the error message

Error, (in Optimization:-LPSolve) problem must have at least one variable

This is not too surprising, because when I examine the arguments being supplied to the LPSolve() command on the first iteration of the for loop these are 0, {}, {} respectively. That is, the LPsolve() command

Optimization[LPSolve](z[Ni], {aconu, cons3, cons4}, binaryvariables = {seq(seq(x[i, j], i = s minus {j}), j = s)})

is actually

Optimization[LPSolve](0, {}, binaryvariables = {})

so I'm not too surpised that this fails miserably.

For a more detailed commentary on why this occurs, see the attached

Solving CVRP
(Fuzzy C-Means Clustering - TSP Problem)

>

restart:

>

gc(): # Force collection
stgc:= kernelopts(gctimes):
st:= time():
FuzzyCMeans:= proc(
     X::Matrix,   #The rows of X are the data points.
     #number of clusters and number of means:
     {clusters::And(posint, satisfies(x-> x > 1)):= 2},
     {fuzziness::satisfies(x-> x > 1):= 2},   #fuzziness factor
     {epsilon::positive:= 1e-4},   #convergence criterion
     {max_iters::posint:= 999},   #max number of iterations
     {Norm::procedure:= (x-> LinearAlgebra:-Norm(x,2))}
)
option `testing`;
description
    "Fuzzy C-Means clustering algorithm.

    "
;
uses LA= LinearAlgebra;
local
     c:= clusters,     #number of clusrters and number of menas
     m:= fuzziness,    #fuzziness factor
     n:= op([1,1], X), #number of points
     d:= op([1,2], X), #dimension
     N:= Norm,
     #U[i,j] is the degree (on a 0..1 scale) to which point i belongs to cluster j.
     U:= Matrix((n,c), datatype= float[8]),
     #U1 is just an update copy of U.
     U1:= LA:-RandomMatrix((n,c), generator= 0..1, datatype= float[8]),
     C:= Matrix((c,d), datatype= float[8]), #The rows of C are the cluster centers.
     e:= 2/(m-1),
     i, j, k, jj, #loop indices
     UM # U^~m
;
     for jj to max_iters while LA:-Norm(U-U1, infinity) >= epsilon do
          U:= copy(U1);
          UM:= U^~m;
          for j to c do
               C[j,..]:= add(UM[i,j].X[i,..], i= 1..n) / add(UM[i,j], i= 1..n)
          end do;
          U1:= Matrix(
               (n,c),
               (i,j)-> 1/add((N(X[i,..]-C[j,..])/N(X[i,..]-C[k,..]))^e, k= 1..c),
               datatype= float[8]
          )
     end do;
     if jj > max_iters then error "Didn't converge" end if;
     userinfo(1, FuzzyCMeans, "Used", jj, "iterations.");
     (C,U1)
end proc:

>

>	#Test data(Clustering data based on the level of demand (low, medium and high)): X:= [[20,20,20],[20,20,20],[15,25,30],[15,15,20],[10,10,15] ] : XM:= Matrix(X, datatype= float[8]):

>	infolevel[FuzzyCMeans]:= 1:

>	(C,U):= FuzzyCMeans( XM, #All options are set to their default values here: clusters= 2, fuzziness= 2, epsilon= 1e-4, max_iters= 999, Norm= (x-> LinearAlgebra:-Norm(x,2)) ):

FuzzyCMeans: Used 2 iterations.

>

AssignToClusters := proc (U::Matrix) options operator, arrow; map(proc (i) local p; member(max(U[i, () .. ()]), U[i, () .. ()], p); p end proc, [`$`(1 .. op([1, 1], U))]) end proc

>

(1)

>

a := [0, op(A)]; NC := 6; objf := Matrix(`$`(NC, 2), shape = symmetric, scan = triangular, `~`[[0, op]]([[1, 3, 2, 1, 2], [3, 2, 6, 7], [3, 5, 6], [2, 2], [4]]))

(2)

>

Error, (in Optimization:-LPSolve) problem must have at least one variable

(3)

>

NULL

>

Download FC_comment.mw

What is the issue with infinity?...

Answered: tomleslie 13876

February 18 2019

0 0

The solution of your final execution group is essentially dependent on the sign of 'a', so the overall answer will be either +infinity or -infinity, which is illustrated in the attached.

So what is tha issue?

>

N := int((x-Q)^m*f(x), x = Q .. infinity); eval(N, [m = 2, f(x) = 1/5000])

int((x-Q)^m*f(x), x = Q .. infinity)

(1)

>

a*(int((x-Q)^m*f(x), x = Q .. infinity))

(2)

>

a*(int(-(x-Q)^m*m*f(x)/(x-Q), x = Q .. infinity))

(3)

>

Val := eval(Dcost, [m = 2, f(x) = 1/5000, co = 25, cs = 15])

(4)

>

`assuming`([simplify(Val)], [0 < a]); `assuming`([simplify(Val)], [0 > a])

>

(5)

>

Download infProb.mw

There are (at least) two issues...

Answered: tomleslie 13876

February 18 2019

0 0

one of which is trivlai - the other is "conceptual"

The first problem is that (in the code you supply) 'x0' is a column vector, whereas 'x1' and 'x2' are row vectors, so the calculation of 'lambda' will fail because you can't add a row vector to a column vector. You can get round this issue by writing the definition of 'x0' as

x0:=Vector[row]([1,1,1,1,1,1,1,1,1,1]);

With 'x0', 'x1', 'x2' now all row vectors, the computation

lambda:=beta__0*x0+beta__1*x1+b__2*x2;

will successfully return a row vector. This however leads to the next problem: exactly what is meant by the expression

B:= x->(1/x!)*lambda^x*exp(-lambda);

This can be defined because the 'body' of the function is not evaluated until it is called. However when you to evaluate this function with any supplied argument, what do you expect from the term exp(-lambda)???

Remember that lambda is a vector; Forget all the complicated 'Statistics' stuff for the moment, what would you expect exp(<1,1>) to return?

Another workaround...

Answered: tomleslie 13876

February 18 2019

0 1

Don't evaluate the the call to BesselJ() before the plot() algorithm gets to work, as in

Download BessPlot.mw

When you write...

Answered: tomleslie 13876

February 17 2019

0 0

F(w)=(10/w)*exp(-7iw)*sin(4*w)

I assume you mean

F(w)=(10/w)*exp(-7*I*w)*sin(4*w)

The default character which Maple uses for hte square root of -1 is uppercase 'I'. You also have to ensure that all multiplication is explicit: so -7*I*w rather than -7Iw

See the atteched

>

  restart;
  F:=(10/w)*exp(-7*I*w)*sin(4*w):
#
# Simple amplitude spectrum, with a few plot
# options
#
  plot( abs(F),
        w = 1..10,
        title="ampl vs frequency",
        titlefont= [times, bold, 16],
        labels=["freq", "ampl" ],
        labelfont= [times, bolditalic, 16]
      );
#
# Of course you can plot the amplitude in decibels
# and the frequency on a logarithmic axis if
# desired/necessary. Shown below with a few more
# plot options
#
  plots:-semilogplot( 20*log10(abs(F)),
                      w=1..100,
                      axes=boxed,
                      gridlines=true
                    );

>

Download spectrum.mw

An important concept you have to underst...

Answered: tomleslie 13876

February 16 2019

1 1

I'm only adding this because I don't think the use of unapply() fully explains your underlying concepts which you really have to get your head around.

Consider the 'toy' example

x:=1;
y:=2;
f:=(x,y)->x+y

So what is f(3,4)???

Think about it, because it is important!

Maple will return 7, because the values used within the function 'f' are local to the function f. The fact that 'x' and 'y' have also been assigned outside the function 'f' is completely irrelevant. Within the function f these assignments do not exist.

Now consider your actual problem

aux := rsolve({y(0) = y0, y(n) = 4*y(n-1)*(1-y(n-1))}, y(n));
solucao := (n,y0)-> aux;
solucao(3,1/2);

Corresponding to the toy example which I gave above, the values 'n' and 'y0' which are used/defined within 'aux' bear no relation to the arguments 'n' and 'y0' passed to the function 'solucao'. Hence they will not be used within the expression assigned to 'aux'. Hence the call solucao(3,1/2) will return an expression containing 'n' and 'y0'

The most common but erroneous next attempt is to try

solucao := (n,y0)-> rsolve({y(0) = y0, y(n) = 4*y(n-1)*(1-y(n-1))}, y(n)); ****
solucao(3,1/2);

After all 'y0' and 'n' are now only defined within the function 'solucao', which is true - so everything ought to be good. Errrrr no!

When you define functions/procedures this way, the 'body' of the function/procedure is not evaluated when it is defined, it is just sitting there waiting for some useful parameters. When you supply thse parameters, they are immediately substituted before anything else happens, so in the above example, the first step performed by the function 'solucao' is to substitute the passed arguments, everywhere they occur within the function, thus creating

rsolve({y(0)=1/2, y(3)=4*y(3-1)*(1-y(3-1))}, y(3))

and I'm sure you can see that rsolve() cannot create a general closed-form recursion relation based on this.

So why does unaspply() work? Firstly the arguments are localised, so the first issue I described above does not occur. Secondly the arguments to unapply() are evaluated, unlike the result of the statement marked **** above. Since the rsolve() command is evaluated, it produces a simple expression containing 'n' and 'y0' which will return a value when the passed arguments are substituted. You should compare the outputs of the two commands

solucao:= (n,y0)-> rsolve({y(0) = y0, y(n) = 4*y(n-1)*(1-y(n-1))}, y(n));
solucao:= unapply(rsolve({y(0) = y0, y(n) = 4*y(n-1)*(1-y(n-1))}, y(n)), n, y0);

to appreciate the difference

Sometimes...

Answered: tomleslie 13876

February 16 2019

0 0

it is more effective to "simplify" on a term-by-term basis.

See the attached

Download simpl.mw

Insufficient precision...

Answered: tomleslie 13876

February 16 2019

2 4

I reran your worksheet with Digits=30 and what I assume to be your "residuals" (ie final two columns of output table) decrease to ~10^-12 or less

I also added the options 'complex' and 'fulldigits' to the fsolve() command - these may be superfluous!

See the attached

>

(1)

>

(2)

>

low_bnd := 0; hig_bnd := 10^10+I*10^10

>

$for b to 20 do Grl_Dst_Frq_Mid_var := fsolve(EQ1 = 0, complex, fulldigits, omega = low_bnd .. hig_bnd_Mid); if `or`(type(Grl_Dst_Frq_Mid_var, complex), type(Grl_Dst_Frq_Mid_var, float)) then Frq_Rate := evalc(Im(Grl_Dst_Frq_Mid_var)/Re(Grl_Dst_Frq_Mid_var)); EQV := evalf(subs(omega = Grl_Dst_Frq_Mid_var, EQ1)); printf("b=%2d ω=%16.8z EI R=%10.8f eqv=%16.8z EI\n", b, Grl_Dst_Frq_Mid_var, Frq_Rate, EQV); hig_bnd_Mid := `if`(b = 20, hig_bnd, Grl_Dst_Frq_Mid_var) else Grl_Dst_Frq_Mid_var := 0.; Frq_Rate := 0.; b := 20 end if end do$

b= 1 ω= 4.99987883E+09 5.00011822E+09I R=1.00004788 eqv= 6.10172162E-12 7.57805552E-12I

b= 2 ω= 2.49982154E+09 2.50024927E+09I R=1.00017110 eqv= 3.05407246E-14 -1.15033788E-13I
b= 3 ω= 1.24959759E+09 1.25035092E+09I R=1.00060285 eqv= 1.03183991E-15 1.65489258E-15I

b= 4 ω= 6.24306563E+08 6.25608647E+08I R=1.00208565 eqv= 4.67983829E-20 -1.61171221E-18I
b= 5 ω= 3.11319803E+08 3.13512118E+08I R=1.00704200 eqv= 1.94136600E-20 -1.56455163E-20I

b= 6 ω= 1.54525508E+08 1.58056231E+08I R=1.02284881 eqv= 4.49552329E-22 3.72286744E-22I

b= 7 ω= 7.61164482E+07 8.13195602E+07I R=1.06835726 eqv= 4.65075796E-24 -1.45859676E-24I
b= 8 ω= 3.61078255E+07 4.29264729E+07I R=1.18884127 eqv= -9.47456282E-26 -2.94877892E-26I

b= 9 ω= 1.58885629E+07 2.32993548E+07I R=1.46642304 eqv= 0.00000000E+00 0.00000000E-01I
b=10 ω= 7.17644439E+06 1.09487978E+07I R=1.52565772 eqv= -1.27217286E-27 -5.79732121E-28I

>

Download moreDigits.mw

The attached works...

Answered: tomleslie 13876

February 16 2019

0 0

although I'm not 100% sure why!!!

So what did I do?

independent variables appeared sometimes as 'bare' names and sometimes as 'backquoted' names. This probably(?) shouldn't be an issue, but I wanted to eliminate anything I didn't "like", so I renamed the independent variables to xi__1 and xi__2. This renaming involves converting the whole worksheet to 1-D input, becuase I have never worked out how to make something as simple as find_and_replace work in 2-D Math. I doubt(?) if this change did anything significant
Your ODEs contain some pretty horrible numerical parameters - I mean 2.751650679*10^(-360)- really!! This is actually too small to be expressed in hardware floats, which are limited to about 10^(-308) roughly. It is also true that many of your numeric parameters are quoted to at least 10 digits which is the default length for Maple's software floats. For both of these reasons I decided to "force" the use of software floats by setting Digits:=30. It is also true that something which is a "numeric singularity" at 10-digit precision, may not be a "numeric singularity" if more precision is used.

Anyhow, with the above changes, the attached works. If I had the time/patience, I might experiment with the setting of Digits to see exactly what level of precision makes the singularity error appear/disappear

>	restart; Digits := 30; kernelopts(version);

(1)

>

fff1 := -(6546.295395*(2.302290809+3.164903854*r^4-5.043972799*r^2-1.539281515*r^6+.5795250834*r^8))*xi__1(r)+((-982.1738168*r-1350.816146*r^5+2152.823658*r^3+656.9824609*r^7-247.3477475*r^9)*cos(1.920758532*r)+(2.379420223*10^(-14)*r+3.272495356*10^(-14)*r^5-5.215443603*10^(-14)*r^3-1.59160968*10^(-14)*r^7+5.992261477*10^(-15)*r^9)*cos(5.762275596*r)+(1.907730247*10^(-99)*r+2.623764526*10^(-99)*r^5-4.181547849*10^(-99)*r^3-1.276093183*10^(-99)*r^7+4.8043714*10^(-100)*r^9)*cos(13.44530972*r)+(-1.043558153*10^(-167)*r-1.435240055*10^(-167)*r^5+2.28737179*10^(-167)*r^3+6.980428435*10^(-168)*r^7-2.628065972*10^(-168)*r^9)*cos(17.28682679*r)+(4.382640639*10^(-253)*r+6.02759067*10^(-253)*r^5-9.606296047*10^(-253)*r^3-2.931576858*10^(-253)*r^7+1.103711249*10^(-253)*r^9)*cos(21.12834385*r)+(-1.437372622*10^(-355)*r-1.976866121*10^(-355)*r^5+3.150572466*10^(-355)*r^3+9.614679048*10^(-356)*r^7-3.619836672*10^(-356)*r^9)*cos(24.96986092*r)+(-2.586708619*10^(-48)*r-3.557585944*10^(-48)*r^5+5.669798376*10^(-48)*r^3+1.730266237*10^(-48)*r^7-6.514290437*10^(-49)*r^9)*cos(9.60379266*r)+1.287492742*10^(-52)*sin(9.60379266*r)+.2444306348*sin(1.920758532*r)-1.973863778*10^(-18)*sin(5.762275596*r)-6.782443883*10^(-104)*sin(13.44530972*r)+2.885635218*10^(-172)*sin(17.28682679*r)-9.915405045*10^(-258)*sin(21.12834385*r)+2.751650679*10^(-360)*sin(24.96986092*r))/r+(34.20619008*(diff(xi__1(r), r)+r*(diff(xi__1(r), r, r))))/r+(226.3864608*(diff(xi__2(r), r)+r*(diff(xi__2(r), r, r))))/r:

>

fff2 := -(11.33151747*(2.302290809+3.164903854*r^4-5.043972799*r^2-1.539281515*r^6+.5795250834*r^8))*xi__2(r)+(-.817717618*r*cos(1.920758532*r)^2+(-27.62404336*r-38.13976873*r^5+60.78413906*r^3+18.54964438*r^7-6.983767491*r^9)*cos(1.920758532*r)-.124054546*r*cos(5.762275596*r)^2+(1.077334797*r+1.487446985*r^5-2.370575056*r^3-.7234341874*r^7+.2723662005*r^9)*cos(5.762275596*r)-0.2342489574e-1*r*cos(13.44530972*r)^2+(0.8525383383e-1*r+.11770766*r^5-.1875931349*r^3-0.5724825577e-1*r^7+0.2155343248e-1*r^9)*cos(13.44530972*r)-0.1420509071e-1*r*cos(17.28682679*r)^2+(-0.4013141059e-1*r-0.5540835201e-1*r^5+0.8830543775e-1*r^3+0.2694838641e-1*r^7-0.1014581525e-1*r^9)*cos(17.28682679*r)-0.9520879006e-2*r*cos(21.12834385*r)^2+(0.2198553924e-1*r+0.3035483874e-1*r^5-0.4837713496e-1*r^3-0.147633686e-1*r^7+0.5558270095e-2*r^9)*cos(21.12834385*r)-0.6821537113e-2*r*cos(24.96986092*r)^2+(-0.1332123502e-1*r-0.1839226852e-1*r^5+0.2931213911e-1*r^3+0.8945257182e-2*r^7-0.3367805606e-2*r^9)*cos(24.96986092*r)-0.4564157749e-1*r*cos(9.60379266*r)^2+(-.2336697743*r-.3226215304*r^5+.5141686128*r^3+.156910093*r^7-0.5907518141e-1*r^9)*cos(9.60379266*r)-.3086246354*r+0.1061426804e-3*sin(9.60379266*r)+0.6274003595e-1*sin(1.920758532*r)-0.8156182761e-3*sin(5.762275596*r)-0.2766135142e-4*sin(13.44530972*r)+0.1012743424e-4*sin(17.28682679*r)-4.539436579*10^(-6)*sin(21.12834385*r)+2.327333679*10^(-6)*sin(24.96986092*r))/r+(15.99132273*(diff(xi__2(r), r)+r*(diff(xi__2(r), r, r))))/r-1.922675548*(diff(xi__1(r), r))*(diff(xi__2(r), r))+(7.650340748*10^(-21)*sin(5.762275596*r)-4.990090145*10^(-55)*sin(9.60379266*r)+2.628753179*10^(-106)*sin(13.44530972*r)-1.118420275*10^(-174)*sin(17.28682679*r)+3.843032538*10^(-260)*sin(21.12834385*r)-1.066490279*10^(-362)*sin(24.96986092*r)-0.9475055746e-3*sin(1.920758532*r))*(diff(xi__1(r), r))+(-2.760876781*10^(-10)*sin(24.96986092*r)-0.1122658745e-2*sin(1.920758532*r)+1.794770960*10^(-6)*sin(5.762275596*r)-8.478174653*10^(-8)*sin(9.60379266*r)+1.129847090*10^(-8)*sin(13.44530972*r)-2.504750041*10^(-9)*sin(17.28682679*r)+7.519222634*10^(-10)*sin(21.12834385*r))*(diff(xi__2(r), r))-0.1385155998e-1*xi__2(r)*xi__1(r)+(0.1716963707e-3*cos(5.762275596*r)-0.3754920005e-4*cos(9.60379266*r)+0.1373099815e-4*cos(13.44530972*r)-6.469641828*10^(-6)*cos(17.28682679*r)+3.546004779*10^(-6)*cos(21.12834385*r)-2.149143420*10^(-6)*cos(24.96986092*r)-0.3977737145e-2*cos(1.920758532*r))*xi__2(r)+(-0.1175680450e-2*cos(1.920758532*r)+2.283732267*10^(-105)*cos(13.44530972*r)-1.249237114*10^(-173)*cos(17.28682679*r)+5.246432433*10^(-259)*cos(21.12834385*r)-1.720670015*10^(-361)*cos(24.96986092*r)-3.096533143*10^(-54)*cos(9.60379266*r)+2.848389467*10^(-20)*cos(5.762275596*r))*xi__1(r)-14.32155064*(diff(xi__2(r), r))^2+(-5.521753562*10^(-9)*sin(24.96986092*r)-0.2245317490e-1*sin(1.920758532*r)+0.3589541920e-4*sin(5.762275596*r)-1.695634931*10^(-6)*sin(9.60379266*r)+2.259694180*10^(-7)*sin(13.44530972*r)-5.009500082*10^(-8)*sin(17.28682679*r)+1.503844527*10^(-8)*sin(21.12834385*r))*(diff(xi__2(r), r))-.1353682171*xi__2(r)^2+(0.3433927414e-2*cos(5.762275596*r)-0.7509840010e-3*cos(9.603792660*r)+0.2746199630e-3*cos(13.44530972*r)-0.1293928366e-3*cos(17.28682679*r)+0.7092009558e-4*cos(21.12834385*r)-0.4298286840e-4*cos(24.96986092*r)-0.7955474290e-1*cos(1.920758532*r))*xi__2(r):