Carl Love

## 26663 Reputation

11 years, 231 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

## Multi-assignment and multi-indexing...

Your can simplify your procedure substantially while correcting the error that you asked about. In the following, I use the optional argument to list the indices that you want to perform this operation on (because it's not clear to me where your index-pattern 1, 2, 4, 5, 7 comes from). The index -1 refers to the last vector position regardless of the actual number of entries.

oneStep_egcd2:= proc(prev::Vector, curr::Vector, J::list(posint):= [1,2,4,5,7])
if curr[-1] <> 0 then
(prev[], curr[J]):= (curr[], prev[J] - curr[J]*iquo(prev[-1], curr[-1]))
end if
end proc
:
oneStep_egcd2(prev, curr);

Maple's multi-assignments happen essentially simultaneously without needing a temporary intermediary variable. For example, compare these:

(a,b):= (1,2):  (a,b):= (b,a):  a, b;  #multi-assignment
2, 1

(a,b):= (1,2):  a:= b: b:= a:  a, b;   #sequential assignments
2, 2

## A better algorithm...

I think that I have a better algorithm for this. I think it has a much better average-case time. It guarantees return of the minimal possible (if there's any y larger than its corresponding x). It does this by stepping sequentially through y, stopping at the first that's greater than its corresponding x. And, it's a much simpler algorithm theoretically.

MinY:= proc(Coeffs::list(rational), m::list(posint))
local c, M:= ilcm(m), g, X:= 0, Y:= 0;
`mod`:= modp; #Use nonnegative remainders only
M/= (g:= igcd((c:= chrem(Coeffs, m)), M)); #chrem is Chinese Remainder Algorithm
c:= g/c mod M;
to M do until (X:= X+c mod M) < (Y+= g);
`if`(X<Y, [X,Y], FAIL)
end proc
:
# -x_coeff/y_coeff for each equation:
Coef:= [-154/69, -13/716, -23/3059, -2295/4522, -6479/5396]:

# Moduli: The 3rd and 4th have their power reduced because the coeffs are reducible.
M:= [7^3, 13^3, 23^2, 17^2, 29^3]
:
seq(CodeTools:-Usage(MinY(Coef[..k], M[..k])), k= 1..nops(M));
memory used=2.80KiB, alloc change=0 bytes,
cpu time=0ns, real time=2.00ms, gc time=0ns

memory used=3.09KiB, alloc change=0 bytes,
cpu time=0ns, real time=0ns, gc time=0ns

memory used=12.67KiB, alloc change=0 bytes,
cpu time=0ns, real time=0ns, gc time=0ns

memory used=16.86KiB, alloc change=0 bytes,
cpu time=0ns, real time=0ns, gc time=0ns

memory used=95.96MiB, alloc change=0 bytes,
+cpu time=2.55s, real time=2.03s, gc time=796.88ms

[7, 49], [691, 819], [15893, 18837], [1103, 26390], [186957261, 190687133]

## Array indices must be integers...

Array indices must be integers. You're trying to use Pi/3 as an index.

## The power series...

A fairly simple power series for the incomplete upper Gamma function is

GAMMA(a,z) = GAMMA(a) - z^a*Sum((-1)^n*z^n/n!/(n+a), n= 0..infinity)

provided that is not a nonpositive integer. Although this series involves a possibly noninteger power of z, that can be factored out of the power series proper, as shown.

## MemoryInUse alternative...

Here's my alternative version of MemoryInUse that runs 68 times faster on my computer:

MemUse:= proc() local r; add(r[3], r= kernelopts('memusage')) end proc:

This measures exactly the same thing as MmaTranslator:-Mma:-MemoryInUse. There is something akin to a "Heisenberg Uncertainity Principle" involved here: It's impossible to measure memory usage without having some effect on it. And of course that's compounded by the fact that the device doing the measuring is the same as the device being measured. So don't expect identical numbers.

It's not incorrect for the memory usage number to sometimes decrease. This is due to "garbage collection"---the recycling of memory addresses whose contents the kernel determines will no longer be needed. Every time that you see a change in the "Memory: " or "Time: " number on the "status bar" (bottom border of your Maple window), a garbage collection has occurred.

Although there is some possibility that "remember tables" (which are the things that the forget command forgets) are causing your memory trouble, I think that that chance is only slight. So you should focus on other causes for now.

I think that you've gotten the mistaken idea that remember tables are primarily used with trig functions. Actually, any piece of Maple code could involve hundreds of remember tables, most of which will be associated with procedures that the user didn't directly call. But, most such remember tables are cleared at every garbage collection (this is called option system).

## No equivalence operator...

The Maple language doesn't have a builtin operator for logical equivalence. I guess that the 2D Input parser translates a variety arrow-type logical symbols into the language's implies operator (which prettyprints as a character like =>).

But, any character that has an HTML code (and thus begins with `&...in ASCII) can be used as an inert infix operator, so what you want prettyprinted can be achieved by

`&iff;`(not (P and Q), ``(not P) or ``(not Q));

The quotes on ``(not P) and ``(not Q) are needed to prevent automatic simplification.

You may want to get rid of some parentheses from the above (at which point the expression will only be good for display, not computation). That can be done like this

use T=Typesetting in
T:-mrow(
T:-Typeset(not (P and Q)),
T:-mo(`&iff;`),
T:-Typeset(``(not P) or ``(not Q))
)
end use;

## Linear regression...

The problem initially seems nonlinear, but it can be easily converted (completely) to a linear-regression problem. By "completely", I mean that the following is the complete original problem, not just some linear approximation of it. The linear form gives you much more control.

 > restart: Digits:= 15: St:= Statistics :
 > Zscore:= proc(V::Vector) uses St=Statistics; local mu:= St:-Mean(V), sigma:= St:-StandardDeviation(V);     mu, sigma, (V -~ mu)/sigma end proc :
 > #Dependent variable: Hvals:= <210.84, 218.19, 219.17, 222.27, 226.35, 228.45, 230.44, 234.65, 239.75>: H__avg, H__sd, H__z:= Zscore(Hvals) :
 > #Independent variable: Qvals:= <1800, 1750, 1730, 1700, 1650, 1635, 1600, 1570, 1500>: Q__avg, Q__sd, Q__z:= Zscore(Qvals) :
 > n:= numelems(Qvals) :
 > #Proposed model, with (apparently) A known beforehand and Q__0, k__f, k__1, and Q__s to be #determined. model:=     subs(         Q= Q__N*Q__sigma + Q__mu,         H__N = (A*(1 - Q/Q__0) - k__f*Q^2 - k__1*(Q - Q__s)^2 - H__mu)/H__sigma     ) :
 > pv:= St:-LinearFit(     [1,Q,Q^2], Q__z, H__z, Q, summarize,     output= [parametervector, residualsumofsquares] );

Summary:
----------------
Model: .52025574e-1-1.0044218*Q-.58528770e-1*Q^2
----------------
Coefficients:
Estimate  Std. Error  t-value   P(>|t|)
Parameter 1    0.0520    0.0414      1.2565    0.2556
Parameter 2   -1.0044    0.0319     -31.5038   0.0000
Parameter 3   -0.0585    0.0325     -1.8025    0.1215
----------------

Solve for the symbolic paarameters in the model by equating coefficients of the powers of Q. We have 3 coefficients and 4 parameters, so I'll just set the value of k__1 (arbitrary choice) to the MathCAD value, 1e-5.

 > Known_params:=     [A, k__1, Q__mu, Q__sigma, H__mu, H__sigma]     =~ [350.156, 1e-5, Q__avg, Q__sd, H__avg, H__sd] ; PV:= solve(     eval(         {seq(pv[1]=~ coeff~(expand(rhs(model)), Q__N, [0,1,2]))},         Known_params     ) );

Compute the sum of the squared residuals. It should be the same as returned by the LinearFit command, allowing for some round off error.

 > add(     eval(         eval(H__N - rhs(model), [PV[2][], Known_params[]])^2,         [H__N= H__z[k], Q__N= Q__z[k]]     ),     k= 1..n );

Compute the sum of the squared residuals using MathCAD's parameter values.

 > add(     eval(         eval(             H__N - rhs(model),             [                 ([Q__0, k__f, Q__s]=~ [6962.32773, 2e-5, 2198.86103])[],                 Known_params[]             ]         )^2,         [H__N= H__z[k], Q__N= Q__z[k]]     ),     k= 1..n );

Look closely at your solve command:

slove([a11*x+a12*y+a13*z=b1,a21*x+a22*y+a23*z=b2,a31*x+a32*y+a33*z=b3[,[x,y,z]);

Surely you can figure out what's wrong with it.

In addition to an `evalf/constant/...` procedure, there are a few more details that you should add to a formal symbolic constant's definition to help it work with the type system and the assume/property/is system.

restart:
`evalf/constant/K`:= ()-> evalf(sqrt(2)*sin(sqrt(17))):
constants,= K;  #Or constants:= constants, K
protect(K);

## InertForm:-Display...

I don't have much familiarity with embedded components, but this may do what you want:

InertForm:-Display(
`%*`(seq(`if`(p[2]=1, p[1], p[1]%^p[2]), p= ifactors(2600)[2])),
inert= false
);

Note that I used the command ifactors instead of ifactor.

Let me know if this helps.

## Your 1st example's rank is 2...

An matrix with a column of zeros (such as your first example) cannot possibly have rank n, regardless of the algebraic field used. The rank of both of your example matrices is 2.

## Multiple initial condition syntax...

Change [x(0)=1, y(0)=1] to [[x(0)=1, y(0)=1]]. The reason is that DEplot allows you to specify multiple sets of initial conditions in a single plot, so the syntax requires them to be subgrouped.

Also, if you haven't loaded package DEtools with with(DEtools), then change DEplot to DEtools:-DEplot. I strongly prefer this package:-member syntax instead of with when one doesn't need the whole package.

I get this plot

## =~...

You almost had it. As you figured, the fundamental command for this is combinat:-randperm. So, you permute the right sides of the equations, and then match them with the left sides. The matching is most easily done with =~ (the elementwise equation-building operator):

([f,g,h,i] =~ combinat:-randperm([f||(1..4)]))(x);
[                  x          2         1]
[f(x) = x, g(x) = 3 , h(x) = x , i(x) = -]
[                                       x]

## List indexing...

It can be done like this:

alpha:= [15, -5/2, 15/11, -15/16, 5/7]:
p:= 3:  #starting position
beta:= alpha[[seq(p..1, -1), p+1..-1]];

I'm not sure if this'll work in your older Maple. If it doesn't, let me know; there'll certainly be something similar that works.

This is probably obvious: The elements of alpha can be anything; they needn't be numbers, real or otherwise.

## ':-p' < ':-c'...

Use ':-p' < ':-c'. The purpose of the :- is to make it refer to the global variable rather than the parameter. The purpose of the single quotes is to make it ignore any values that may be assigned to the globals. Unevaluation techniques (such as single quotes) usually don't work with procedure parameters, which is why you need to switch to the globals.

 First 8 9 10 11 12 13 14 Last Page 10 of 384
﻿