Carl Love

Carl Love

28020 Reputation

25 Badges

12 years, 301 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are replies submitted by Carl Love

@vsnyder You wrote:

  • Another problem (that isn't mentioned in the Maple Advanced Programming Guide) is that sending statements to EvalMapleStatement isn't quite the same as typing them into the Maple textual interrace. In the Maple textual interface, you can lay out a "proc" definition on several lines, as you would in C or Fortran or Haskell. To submit it to EvalMapleStatement, it needs to be combined into one long stream of text, from "proc" to "end proc".

That is completely false. You either haven't understood anything that I've said in this thread regarding line breaks, white space, etc., or you've chosen to ignore it.

@MapleEnthusiast You have a system of 4 equations in 4 unknowns, so I don't know why you think that its solution has a free variable that you can arbitrarily set to 1. 

It is much, much, much easier to find numeric solutions of sysP with fsolve than it is to work with the huge symbolic solution I generated in the Answer. For example:

fsolve(sysP, {vz[]}=~ 0.1..infinity);
        {z[2, 1] = 1.204997696, z[2, 2] = 1.168035195,  
         z[4, 1] = 3.038404729, z[4, 2] = 3.170903246}

Remember that these need to be squared to get the corresponding y values.


 

@elonzyy You wrote:

  • to_list := a -> convert(a, list);
    

Better (for generic a): 

to_list:= `convert/list`:

If you're expecting a to be a certain type (Vector, table, set, etc.), then there are better ways.

  • map_i := (f, a) -> zip(f, to_list(a), [`$`(1..nops(to_list(a)))]);

Yes, that's bloated due to unnecessary quotes and parentheses and inefficient due to the repetition of to_list. Better:

map_i:= (f,a)-> (L-> f~(L, [$1..nops(L)]))(to_list(a)):

  • flat_map_i := (f, a) -> ListTools:-Flatten(map_i(f, a), 1);
    flat_map := (f, a) -> ListTools:-Flatten(map(f, a), 1);

Those two aren't even worth typing or giving a name to. You use an f that creates a sublist from a sequence, and then you turn that sublist back into a sequence. It'd be better to simply not create any list with in the first place, for example,

map_i(1@@0, [3,5,7]);

(1@@0 is the multi-argument identity function. An equivalent form for it is ()-> args.)

 

@Zeineb The reason that your instructor assigned you this exercise is specifically to teach you some of the same things that I've been telling you. This function was specifically designed to cause trouble (an erratic sequence) for Newton's method when the initial value is several "humps" away from the true solution. Look at a plot of f and another of its derivative on the interval x= 0..22. Note the wild oscillations of the derivative. If the derivative is near zero, the tangent line is nearly horizontal, and thus the point where the tangent intersects the x-axis (which is, by definition, the next x) will be far away from the current x.

Here's a procedure that I think will help you see the "bare bones" of what's happening. You can use it as an alternative to Newton in Student:-NumericalAnalysis. You should especially experiment with different values of digits and different initial values.

Newton:= proc(
    f::algebraic, 
    X0::(name=realcons),
    { 
        digits::posint:= Digits,
        tolerance::positive:= .5*10^(2-Digits),
        maxiters::posint:= 20,
        criterion::identical(relative, absolute, residual):= ':-relative'
    }
)
local 
    X:= Array(1..2, [rhs(X0)+1, rhs(X0)], datatype= sfloat),
    x:= lhs(X0),
    F:= subs(
        _f= unapply(f,x), 
        proc(x) option remember; evalf(_f(x)) end proc
    ),
    N:= subs(_d= unapply(diff(f,x), x), x-> x-F(x)/evalf(_d(x))),
    Crit:= table([
        ':-absolute'= ((a,b,t)-> abs(a-b) < t),
        ':-relative'= ((a,b,t)-> abs(1-b/a) < t),
        ':-residual'= ((a,b,t)-> abs(F(a)) < t)
    ]),
    k
;
    Digits:= digits;
    for k from 3 to 2+maxiters do
        if Crit[criterion](X[k-1], X[k-2], tolerance) then break fi;
        X(k):= N(X[k-1])
    od;
    if k > 2+maxiters then print("Did not converge") fi;
    seq(X[2..])
end proc
:   
f:= sin(2*x)+3*x+2*cos(x)^2-61:
Newton(f, x=0, digits= 15, tolerance= .5e-13, criterion= relative);

 

@vsnyder I meant a transcript of the Fortran calls to EvalMapleExpression, not a transcript of how the results appear  within Maple. But at this point I think that it's clear to both of us what you were doing wrong: The procedure definition, or any expression, must appear as a single call to EvalMapleExpression. It's make no difference whether the argument to that call contains line breaks. Analogously, if you were entering code into Maple by hand, regardless of interface, no single expression (even a thousand-line module) could be split over multiple execution groups (command prompts); nor would you be allowed to press Return until it was completely entered. 

The file attachment tool in MaplePrimes is the green uparrow rather than the more-usual paperclip icon on a lot of sites. At least the bright color makes it easier to see.

@Rookieplayer Yes, I noticed many SLATEC routines related to piecewise polynomials (just use your browser's in-page search feature (usually Ctrl-F) on the page https://gams.nist.gov/cgi-bin/serve.cgi/PackageModules/SLATEC), but nothing for non-polynomials. Of course, your functions are piecewise-analytic, so they could be approximated by piecewise polynomials; however, I think that the approach that I described in my previous Reply will be far more accurate, far more computationally efficient, and much easier to implement.

By the way, in case you don't already know this: You should check out Maple's command CodeGeneration:-Fortran, which I used to generate HVS above.

@Rookieplayer I can't find piecewise in SLATEC either. Before sending your code to Fortran, convert all the piecewises to Heaviside (this could be done all at once with the subsindets command). I can't find Heaviside in SLATEC either, but it'd be trivial to implement. Substitute HVS (or another Fortran-acceptable name) for Heaviside (this could be done all at once with subs). Then use Fortran code such as this for HVS:

doubleprecision function HVS (x)
    doubleprecision x
    if (x .le. 0.0D0) then
         HVS = 0.0D0
         return
    else
         HVS = 0.1D1
         return
    end if
end

 

My latex(...results are quite different from yours:

restart:
Physics:-Version();
         The "Physics Updates" package is not installed

interface(version);
Standard Worksheet Interface, Maple 2021.0, Windows 10, March 5 
   2021 Build ID 1523359

interface(typesetting);
                            extended
interface(prettyprint);
                               2
latex:-Settings();
[cacheresults = true, commabetweentensorindices = false, 
  invisibletimes = " ", leavespaceafterfunctionname = false, 
  linelength = 66, powersoftrigonometricfunctions = mixed, 
  spaceaftersqrt = true, usecolor = true, 
  usedisplaystyleinput = true, useimaginaryunit = I, 
  useinputlineprompt = true, userestrictedtypesetting = false, 
  usespecialfunctionrules = true, 
  usetypesettingcurrentsettings = false]

sol:= u(r,t) = 
    invlaplace(
        BesselJ(0,10*(-s)^(1/2)*r)/BesselJ(0,20*(-s)^(1/2))*s/(s^2+1),
        s, t
     ) 
     - invlaplace(
           BesselJ(0,10*(-s)^(1/2)*r)/BesselJ(0,20*(-s)^(1/2))/s,
           s, t
     )
     - cos(t) + 1
:
latex(sol);
u \! \left(r , t\right) = 
\mathcal{L}^{-1}\! \left(\frac{J_{0}\! \left(10 \sqrt{-s}\, r \right)
 s}{J_{0}\! \left(20 \sqrt{-s}\right) \left(s^{2}+1\right)}, s , t\right)
-\mathcal{L}^{-1}\! \left(\frac{J_{0}\! \left(10 \sqrt{-s}\, r \right)}
{J_{0}\! \left(20 \sqrt{-s}\right) s}, s , t\right)-\cos \! \left(t \right)+1

 

I can confirm that I also get the error, exactly like you described and your picture shows. I have no idea what causes it though.

@vsnyder You wrote:

  • I found that if I want to define a proc, it has to be sent to Maple as one line in EvalMapleExpression.

Would you please post an exact transcription of at least two examples (one positive and one negative) of some code and the resulting error message (in the positive case) that provides evidence of the truth of that statement? I'm not denying that it's true; I just need to see it to understand what's happening.  

  • I tried to break it into several lines with backslash at the end, but that resulted in "unexpected end of statement."

And transcription of that too, please. I'd guess that the string and line-breaking rules of Fortran play a role in this also. So, could you summarize those rules?

Consider this example (and perhaps this explains your problem):

"This is a\
multi-line Maple \
string";
               "This is amulti-line Maple string"

Note the transcription error that is caused by ending with just backslash rather than space-backslash. If that had been functional Maple code rather than just a text string, the running together of two words would likely cause a syntax error. A line break (aka a newline or \n) is considered to be white space syntactically (i.e., functionally equivalent to a space). If that character is escaped with backslash, its syntactic function needs to be replaced by a space.

There is a bug in MaplePrimes such that it doesn't properly upload *.maple files. A workaround is to convert it to a *.zip file and upload that.

@pauldaas Here's an improvement that I should've mentioned earlier:

f:= unapply(piecewise(for i to n+1 do x <= a[i], y[i](x) od), x);

Like f:= x-> ..., this also creates an arrow expression procedure. The difference is that unapply​​​ builds the piecewise expression first (and only once), then turns it into a procedure. The f:= x-> ... form rebuilds the piecewise for each value of x.​​​​

@pauldaas Since you had the right idea about embedded for loops, I should show precisely where your syntax was wrong. But first, note that embedded for loops can only be used in 1D input (plaintext input). So, here's your code, line by line:

  1. f(x)=piecewise(

  2. x<a[1],y(x)[1],

  3. for i from 2 to n+1,

  4. a[i-1] <=  x  <=  a[i], y(x)[i];

  5. od:);

Line 1: To define a function (formally called a procedure or arrow expression), line 1 should be

f:= x-> piecewise(

Line 2:

  1. Because of what I say about line 4 below, it's not necessary to separate out the first branch as special. 
  2. An indexed function is y[1](x), not y(x)[1].

Line 3: 

  1. If line 2 is included (although it's not needed), then the embedded for loop (in its entirety) must be enclosed in parentheses.
  2. Replace the comma at the end with the word do.

Line 4:

  1. Chained inequalities such as a < b < c are not allowed in 1D input. If such is needed (although it's not needed in this case), you could use And(a < b, b < c).
  2. Although the And(...form inequality is allowed for the conditions in piecewise, it's not needed, and it's preferable not to use it. The reason is that the piecewise conditions are processed left to right. If condition k is being evaluated, then it's necessarily true that conditions 1 thru k-1 have already been determined to be false.
  3. The semicolon at the end is optional, and in the case of an embedded for loop, I think that it's undesirable (because the values returned to the sequence by the loop necessarily come from the last statement of the loop). 

Line 5:

  1. An embedded for loop is an expression, not a statement, so it's not allowed for there to be a colon or a semicolon immediately after od.
  2. However, the semicolon after the closing parenthesis is fine because it terminates the statement that began with f:=.

 

@vsnyder Yes, I should've said (and indeed, I was just editing my previous Reply to say it) that the penalty that I was referring to was additive, not multiplicative. So, if it cost an extra second, then that's approximately what it'll always cost. If one were loading LinearAlgebra to do a trivial-time operation (such as solving a linear system in hardware floats), then 1 second would be a huge penalty.

The primary reason that I don't use with is that I think that it reduces code readability. Instead I assign the package prefix to a short variable name (such as LA:= LinearAlgebra) and then just use that. 

@vsnyder You are correct that you don't want to be opening a new Maple session repeatedly, and I don't know how to help you with StartMaple. However, it's not "necessary" to load packages with with; and, indeed, doing so incurs a significant performance penalty (for the loading itself; no penalty for using the package's commands). Just use the so-called "long forms" of the command names, such as LinearAlgebra:-LinearSolve.

There was no confusion in what you wrote; however, to avoid possible confusion in the future: Things such as LinearAlgebra are called packages, not libraries. The distinction is important because Maple also has libraries, which are the level of "compiled" code containerization one level larger than packages.

First 121 122 123 124 125 126 127 Last Page 123 of 708