Does it help, if you also supply the optional argument 'WRAPPER' in your define_external call?

You may need a C compiler, for that to work. See the ?WRAPPER help-page .

acer

See the help-page ?if .

In programming, you can use the if..then construction as follows,

if a>b and b>c then
  4;
end if;

There is also a function, `if`.

> `if`(a>b and b>c, 4, NULL);
                         if(b < a and c < b, 4, NULL)
 
> eval(%, {a=7,b=3,c=-1});
                                       4
> eval(%%, {a=2,b=3,c=-1});
                                     NULL

acer

It would not be good if Maple actually assigned to Kvr and Kz15 as part of the solving process. It would be too heavy handed, and prevent other use of those names.

You can use the 2-argument form of the eval() command to evaluate other expressions using the result from fsolve (or solve). Think of it as "substituting in" the values of the solution from fsolve.

> Kz24 := (2.4/1.5)*Kz15;
                           Kz24 := 1.600000000 Kz15
 
> eval(Kz24, soln_z15);
                                  323614.6206

nb. You can get a similar effect using subs(), but I don't recommend it here. This is mathematical evaluation, not substitution into a data structure.

acer

I expect that there are lots of ways to do it.

LtoP := L -> add( L[i]*x^(i-1), i = 1..nops(L) );

LtoP([1, 2, 3]);

Or, for Horner's form,

LtoP := L -> foldr((a,b)->x*(a+b),0,op(L[2..-1]))+L[1];

LtoP([1, 2, 3]);

acer

Could you upload the data to the mapleprimes File Manager? Could you post or upload the precise code that you tried, which gave rise to the result that you quoted?

acer

It sounds to me as if you are asking about facilities for doing abstract linear algebra in Maple. If that's right, then I'm afraid that the answer may be: not much.

acer

If you are trying hard to work out some tricky math and maple programming technique, then the last thing you need is additional hardship due to using 2D Math input and/or Document mode. I would suggest that, until you can nail down the programming and math difficulties, you use 1D input in a Worksheet.

acer

You can do this using Maple's Re() or Im() routines, but it helps if you can also inform Maple that the other quantities H1(k), etc, are purely real.

There are a few ways to do that. One way uses the fact that Maple's evalc assumes that names are generally real.

C(k) := H1(k)/(H1(k) + I*H0(k));

RC:=evalc(Re(C(k)));

IC:=evalc(Im(C(k)));

simplify( RC+I*IC - C(k) );

A second way is to utilize assumptions on those other names via Maple's `assuming` facility.

C := H1/(H1 + I*H0);

Re(C) assuming real;

Im(C) assuming real;

acer

There were two coefficients of air resistance given: one for before deployment and one for after.

If you could alter the first, then the (terminal) velocity prior to deployment could be lower. If you could alter the latter, then the rate of change in velocity following deployment would be greater, and also the steady descent rate would be less.

Skydivers keep the former coefficient near optimal by falling spead-eagle instead of with knees tucked up in cannonball position. After deployment, the size of the parachute is the big factor in air resistance. Maybe you could show graphs of k1 and/or k2 versus impact velocity (if, say, the rip cord is pulled after terminal velocity were attained).

acer

There are a few misconceptions about the contents of the ?EffNumLA help-page apparent in your posting.

• B:=evalf(A):

  The above does not produce a Matrix with a datatype directly suitable for external calling to the fast compiled routines for numerical linear algebra. Instead all it does it produce a Matrix with floating-point entries. See example code below.
• DiagonalMatrix(E) won't work, as DiagonalMatrix expects a Vector or list and you have E as a scalar. I tried to guess as what you might have wanted in the code below, and you can replace in a suitable manner.
• MatrixAdd(B+M) is likely not what you meant, if you intended to add Matrices B and M. You may instead have meant MatrixAdd(B, M).
• C:=MatrixAdd(B, M), "Should overwrite old matrix". No, it will reassign a new Matrix to name C. The old assigned value (a Matrix, say) of C will then become garbage which takes Maple longer to manage and collect. If you really want to re-use the same memory space of the old C then use one of the float datatypes and the inplace option for MatrixAdd. See code below. (Strictly speaking it will re-use the same memory space for the data only in the hardware float[8] or complex[8] datatype case, and in the software sfloat or complex(sfloat) datatype case will re-use only the container Matrix.)
• Numerically, it is better to do a LinearSolve than to multiply by an explicitly formed Matrix inverse. A little transposition footwork can allow this even when the eqautions look "back-to-front" from what you might expect.
• As you can see in the comments, there are a set of progressively better ways to get DD.
• The datatypes suitable for direct external calling are float[8], sfloat, complex[8], and complex(sfloat). Depending on Digits, the datatypes float and complex(float) will be converted to one of those accordingly.

I did not do anything about replacing the LinearAlgebra:- routines with their more efficient programming forms whose names begin with LinearAlgebra:-LA_Main . That too should be described in the ?EffNumLA help-page.

Below is a modified version of your basic code, which might help illustrate some techniques. You might want to force it through just one iteration, to check that my methods of computing `out` (or DD, as I called the final result) are equivalent.

You didn't say whether your data was symmetric. If it were, then you could get by with 'datatype'='float' throughout, and Transpose instead of HermitianTranspose.

You could likely also make LinearSolve act in-place on its second (Matrix) argument, and the Hermitian Transpose call on that argument as well. That seemed to be the best way to get the final result (`out`, or `DD` when acting in-place) , and is what I've done below (with alternatives commented out).

Digits:=trunc(evalhf(Digits)):
#Digits:=32:
N:=16:
with(LinearAlgebra):
 
A:=RandomMatrix(N,generator=-0.1..0.1,'outputoptions'=['datatype'='float']):
B:=Matrix(op(1,A),(i,j)->evalf(A[i,j]),'datatype'='float'):
 
(st,ba,bu):=time(),kernelopts(bytesalloc),kernelopts(bytesused):
for E from 0.01 by 0.0001 to 1.0 do
  C:=MatrixAdd(B,
               DiagonalMatrix(RandomVector(N,generator=-0.1..0.1,
                                           outputoptions=['datatype'='float']),
                                           outputoptions=['datatype'='float']),
               'inplace'=true):
  (l,L):=Eigenvectors(C):
  #Linv:=MatrixInverse(L):
  #DD:=L.MatrixExponential(DiagonalMatrix(l,outputoptions=['datatype'='complex'('float')])).Linv;
  # But why not this, since the MatrixExponential of a Diagonal Matrix
  # is easy.
  #DD:=L.DiagonalMatrix(map(t->exp(t),l),outputoptions=['datatype'='complex'('float')]).Linv;
  # And why not this, further eliminating a multiplication step.
  #DD := HermitianTranspose(LinearSolve(HermitianTranspose(L),
  #                                     HermitianTranspose(Matrix(N,N,(i,j)->exp(l[j])*L[i,j],
  #                                                               'datatype'='complex'('float')))));
  #out:=Eigenvalues(DD):
  DD := Matrix(N,N,(i,j)->exp(l[j])*L[i,j],'datatype'='complex'('float')):
  HermitianTranspose(DD,inplace=true): # Or just form DD with reversed indices.
  LinearSolve(HermitianTranspose(L),DD,inplace=true); # result is in DD!
  HermitianTranspose(DD,inplace=true):
end do:
time()-st,kernelopts(bytesalloc)-ba,kernelopts(bytesused)-bu;

On my machine, in Maple 11, it took 190sec and Maple claimed to have allocated about 10-14 million words, with Digits=14 and N=16, to do all 9900 loop iterations.

With Digits=32 and N=16, and running only 200 times through the loop, then Maple claimed to allocate 80 million words of memory and it was much slower. Tracking memory through the OS (top, with Linux) showed that Maple's mserver actually consumed 125MB of resident memory. You could try the above code, with Digits > evalhf(Digits) , and see if it leaks or runs out of memory while running through all the 9900 loop iterations.

acer

Using normal worked for me, in Maple 11.

> expr := (r*(r+(1/2)*s))*(sin(t)/(r+(1/\
> 2)*s)+(diff(y(x), x))/(r^2+(1/2)*r*s));
                                     /            d         \
                                     |            -- y(x)   |
                                     |sin(t)      dx        |
                 expr := r (r + s/2) |------- + ------------|
                                     |r + s/2    2          |
                                     \          r  + 1/2 r s/

> normal(expr);
                                        /d      \
                             sin(t) r + |-- y(x)|
                                        \dx     /

The simplify command also did the job. It neither of those work for you, you might want to check that you're using Maple's 2D Math's implicit multiplication (ie. no `*`) correctly.

acer

We might be able to show how it can be done, but more detail would likely be necessary for that.

For example, what is Y1? Is it the water liquid mole fraction at stage 1? And what is T1, is it the temperature at stage 1?

The key bit might be to know the full relationship between T1 and {X1,Y1,X0,Y0} with all the governing equations laid out.

As an aside, it may be unlikely that there's a chemist reading here who might have worked on this particular problem. But if the mathematical problem can be abstracted away from the discipline (here, chemistry) then the chances are quite a bit higher that the mathematical Maple programmers here might be able to show the way. So it might be better to try to remove as much of the chemistry from the description as possible, and try to couch the problem more generically. Eg. additional equations must hold -- always, perhaps, and hence at each time step in particular -- and they form an implicit relationship between some of the dependent variables or their derivatives, etc.

acer

I guess that you are already using the ':-output'=':-solutionmodule' option to NonlinearFit?

If you have residualmeansquare assigned to something, then have you tried accessing it as eq_fit:-Results("residualmeansquare") ?

Otherwise, could you post the example? It can be difficult to diagnose the problem without a complete piece of code to reproduce the behaviour.

acer

Check out this routine, for analysis of asymptotes of rational (polynomial) functions like your example.

Student[Precalculus][RationalFunctionTutor]();

To answer your other questions, where we are talking about functions y(x) with the x-axis being the horizontal axis and the y-axis being the vertical axis,

• No. The function sin(x)/x is undefined as x=0, and it does not have a vertical asymptote at that point. A function may be undefined at a point and still not have an infinite limiting value at that point. As the sin(x)/x example shows, the graph needn't be approaching a vertical slope at that point either.
• Yes, a function's graph can cross its horizontal asymptotes. The example of (x-2)/(x^2-5x) does this, by cross its horizontal asymptote (which occurs at y=0, as x tends to either -infinity or +infinity) at the point x=2. No, a function's graph may not cross its vertical asymptotes (else those would not actually be asymptotes, or it would be a graph of a a relation rather than of a function).
• A function's graph can have at most two horizontal asymptotes. Remember that horizontal asymptotes occur according to behaviour as x tends to -infinity or +infinity (which only happens once, for each of those).

acer

What sort of a fit do you want? You mentioned "least squares", but even so there is some choice.

Do you want usual linear regression (minimizing the sum of the squares of distance with one particular component such as z be taken as dependent)?

Or do you want to find a line which minimizes the sum of the Euclidean distance between the data points and that line. That is the othogonal distance regression, a four-parameter problem sometimes solved using singular value decomposition, but which can also be solved by generating an objective from all the data points'  point-to-line distance expressions.

There are some other choices, but those are the common ones. The first might be possible using Statistics[LinearFit] and some footwork. The second might be more "usual".

acer