acer

17418 Reputation

29 Badges

14 years, 138 days

On google groups. (comp.soft-sys.math.maple and sci.math.symbolic)

On stackoverflow.

On math.stackexchange.com.

MaplePrimes Activity


These are answers submitted by acer

Here are two ways. One is a CodeGeneration:-Matlab call on a list of a sequence of equations. The other is a loop of CodeGeneration:-Matlab calls on a list with a single equation.

restart

with(CodeGeneration); with(LinearAlgebra)

A := RandomMatrix(3, 3)

Matrix(%id = 18446884456322477342)

Matlab([seq(seq(KNL[ii, jj] = A[ii, jj], jj = 1 .. 3), ii = 1 .. 3)])

KNL(1,1) = 27;
KNL(1,2) = 99;
KNL(1,3) = 92;
KNL(2,1) = 8;
KNL(2,2) = 29;
KNL(2,3) = -31;
KNL(3,1) = 69;
KNL(3,2) = 44;
KNL(3,3) = 67;

for ii to 3 do for jj to 3 do Matlab([KNL[ii, jj] = A[ii, jj]]) end do end do

KNL(1,1) = 27;

KNL(1,2) = 99;
KNL(1,3) = 92;
KNL(2,1) = 8;
KNL(2,2) = 29;
KNL(2,3) = -31;
KNL(3,1) = 69;
KNL(3,2) = 44;
KNL(3,3) = 67;

 

Download Matlab_-_Resultname_ac.mw

The _d01amc method of evalf(Int...)) can do it, although it needs to be forced to not use evalhf mode (which I do with a trick, by using a list within a procedure). And I split the real and imaginary components , since roundoff error produces small but nonzero imaginary components (artefacts).

I did not try converting the integrand to something more tractable.

The result agrees well with that using a chopped moderate finite interval. (The curve goes to zero quickly, as you saw.)

restart

istar := y*sinh(Pi*y)*exp(-alpha*y^2)*LegendreP(-1/2+I*y, 1+u)

y*sinh(Pi*y)*exp(-alpha*y^2)*LegendreP(-1/2+I*y, 1+u)

alpha := 1; u := 2

2

plot(([Re, Im])(istar), y = 0 .. 40, size = [400, 200])

Hre := proc (Y) []; eval(Re(istar), y = Y) end proc

evalf(Int(Hre, 0 .. infinity, method = _d01amc))

-1.707207851

Him := proc (Y) []; eval(Im(istar), y = Y) end proc

evalf(Int(Him, 0 .. infinity, method = _d01amc))

-0.3469635600e-28

evalf(Int(istar, y = 0 .. 100))

-1.707207851-0.1848178275e-27*I

``

Download Istar_ac.mw

In the lines where you seem to be assigning to a__01 and a__10 you have extra spaces in those names, and you are actually assigning to `a__01 ` and `a__10 ` with 1 extra space at the end of the name.

Also, in the expression you are assigning on the very first line you have an instance of `u__2  ` with 2 extra spaces at the end of that name. That is not the same as the u__2 you use elsewhere.

Below I have removed these extra spaces.

a__10 := (1-u__2)*(u__2-m)-q*(u__2+c)+u__2*(1+m-2*u__2)

(1-u__2)*(u__2-m)-q*(u__2+c)+u__2*(1+m-2*u__2)

(1)

a__01 := -u__2*q

-u__2*q

(2)

a__20 := 1/2*(-6*u__2+2*m+2)

-3*u__2+m+1

(3)

a__11 := -q

-q

(4)

a__02 := 0

0

(5)

b__10 := s

s

(6)

b__01 := -s

-s

(7)

b__20 := -s/(u__2+c)

-s/(u__2+c)

(8)

b__11 := 2*s/(u__2+c)

2*s/(u__2+c)

(9)

b__02 := -s/(u__2+c)

-s/(u__2+c)

(10)

b__03 := 0

0

(11)

a__30 := -1

-1

(12)

a__21 := 0

0

(13)

b__12 := s/(u__2+c)^2

s/(u__2+c)^2

(14)

a__12 := 0

0

(15)

b__21 := (-2*(u__2+c)^2-s)/(2*(u__2+c)^3)

(1/2)*(-2*(u__2+c)^2-s)/(u__2+c)^3

(16)

l__1 := ((-3*Pi)*(1/(2*a__01*(q*s*u__2-s^2))))*(a__10*b__10*(a__02*b__11+a__11^2+a__11*b__02)+a__10*a__01*(a__11*b__02+a__20*b__11+b__11^2)+b__10^2*(a__02*a__11+2*a__02*b__02)-2*a__10*b__10*(-a__02*a__20+b__02^2)-2*a__10*a__01*(a__20^2-b__02*b__20)-a__01^2*(2*a__20*b__20+b__11*b__20)+(a__01*b__10-2*a__10^2)*(-a__11*a__20+b__02*b__11)-(a__01*b__10+a__10^2)*(3*(-a__01*a__30+b__03*b__10)+2*a__10*(a__21+b__12)+b__10*a__12-a__01*b__21))

(3/2)*Pi*(((1-u__2)*(u__2-m)-q*(u__2+c)+u__2*(1+m-2*u__2))*s*(q^2+q*s/(u__2+c))-((1-u__2)*(u__2-m)-q*(u__2+c)+u__2*(1+m-2*u__2))*u__2*q*(q*s/(u__2+c)+2*(-3*u__2+m+1)*s/(u__2+c)+4*s^2/(u__2+c)^2)-2*((1-u__2)*(u__2-m)-q*(u__2+c)+u__2*(1+m-2*u__2))*s^3/(u__2+c)^2+2*((1-u__2)*(u__2-m)-q*(u__2+c)+u__2*(1+m-2*u__2))*u__2*q*((-3*u__2+m+1)^2-s^2/(u__2+c)^2)-u__2^2*q^2*(-2*(-3*u__2+m+1)*s/(u__2+c)-2*s^2/(u__2+c)^2)+(-q*s*u__2-2*((1-u__2)*(u__2-m)-q*(u__2+c)+u__2*(1+m-2*u__2))^2)*(q*(-3*u__2+m+1)-2*s^2/(u__2+c)^2)-(-q*s*u__2+((1-u__2)*(u__2-m)-q*(u__2+c)+u__2*(1+m-2*u__2))^2)*(-3*u__2*q+2*((1-u__2)*(u__2-m)-q*(u__2+c)+u__2*(1+m-2*u__2))*s/(u__2+c)^2+(1/2)*u__2*q*(-2*(u__2+c)^2-s)/(u__2+c)^3))/(u__2*q*(q*s*u__2-s^2))

(17)

map(convert, indets(l__1, name), string)

{"Pi", "c", "m", "q", "s", "u__2"}

(18)

 

Download coeff_ac.mw

Another way,

restart;

 

Multiple angle reduction formula (See here) for integer N.

 

eq:=cos(N*t)=Sum((-1)^k*binomial(N,2*k)*cos(t)^(N-2*k)*sin(t)^(2*k),k=0..N/2);

cos(N*t) = Sum((-1)^k*binomial(N, 2*k)*cos(t)^(N-2*k)*sin(t)^(2*k), k = 0 .. (1/2)*N)

value(eval(eq,[N=6]));

cos(6*t) = cos(t)^6-15*cos(t)^4*sin(t)^2+15*cos(t)^2*sin(t)^4-sin(t)^6


An undocumented implementation of that.

eval("full angle reduction",trigsubs(cos(6*t),annotate));

cos(t)^6-15*cos(t)^4*sin(t)^2+15*cos(t)^2*sin(t)^4-sin(t)^6

 

Download multiple_angle_reduction.mw

To assign to x,

   x := log[2](6);

Note the distinction between the exact symbolic result and the float-approximation.

log[2](6);
                                                                                        
                   ln(6)
                   -----
                   ln(2)

log[2](6.0);                                                                                      

                 2.584962501

What version of Maple are you using, and which platform (Operating System)?

restart;
kernelopts(version);

   Maple 2019.1, X86 64 LINUX, May 21 2019, Build ID 1399874

Optimization:-LPSolve(4*x1+5*x2+5*x3+2*x4,
                      {33*x1 + 49*x2 + 51*x3 + 22*x4 <= 120},
                      maximize, assume ={nonnegint},
                      depthlimit=10); 

          [13, [x1 = 2, x2 = 1, x3 = 0, x4 = 0]]

Optimization:-Maximize(4*x1+5*x2+5*x3+2*x4,
                       {33*x1 + 49*x2 + 51*x3 + 22*x4 <= 120},
                       assume=nonnegint, depthlimit=10);

          [13, [x1 = 2, x2 = 1, x3 = 0, x4 = 0]]

I get the same results using Maple 16.02 for 64bit Linux.

Here is one way to do that.
replaceid_example.mw

Here is another way, using readstat for the original query.
readstat_example.mw

I await the followup, "What I actually want to do is..."

If you want to have the return value be a procedure (or operator) then there's no need to demand that the name of its formal parameter be supplied to LinModel.

Note that by using CurveFitting:-LeastSquares rather than Statistics:-Fit the fitting model can be omitted from the call (as the former computes a linear model in that case).

restart;

Mat := module()
    description "My Package";
    option package;
    export LinModel;
    LinModel := proc(xliste::list,yliste::list,$)
        local x;
        return unapply(CurveFitting:-LeastSquares(xliste,yliste,x),x);
        end proc; # LinModel
    end module: # Mat

L1 := [1,2,4,7];

[1, 2, 4, 7]

L1a := evalf(L1);

[1., 2., 4., 7.]

L2 := [2,4,6,8];

[2, 4, 6, 8]

Mat:-LinModel(L1,L2);

proc (x) options operator, arrow; 5/3+(20/21)*x end proc

func := Mat:-LinModel(L1a,L2);

proc (x) options operator, arrow; HFloat(1.6666666666666676)+HFloat(0.9523809523809522)*x end proc

plots:-display(
  plots:-pointplot(<<L1>|<L2>>,
                   symbol=solidcircle,symbolsize=15),
  plot(func, 0..L1[-1]),
  view=[0..L1[-1],0..func(L1[-1])]
);

 

 

Download CF.mw

Of course, you could still supply the name of the formal parameter to be used for the returned procedure (if you want, but there is no significant distinction -- it's just a dummy name).

I have also changed the type specification of parameter x in the definition of procedure LinModel from ::algebraic to ::name which is tighter (and more correct, IMO).

restart;

Mat := module()
    description "My Package";
    option package;
    export LinModel;
    LinModel := proc(xliste::list,yliste::list,x::name,$)
        return unapply(CurveFitting:-LeastSquares(xliste,yliste,x),x);
        end proc; # LinModel
    end module: # Mat

L1 := [1,2,4,7];

[1, 2, 4, 7]

L1a := evalf(L1);

[1., 2., 4., 7.]

L2 := [2,4,6,8];

[2, 4, 6, 8]

Mat:-LinModel(L1,L2,x);

proc (x) options operator, arrow; 5/3+(20/21)*x end proc

Mat:-LinModel(L1,L2,t);

proc (t) options operator, arrow; 5/3+(20/21)*t end proc

 

Download CF2.mw

It may have occurred to you that this a long-winded way to obtain functionality that is almost directly available from a stock Maple command. But perhaps you are trying to learn the ins and outs of packages, and this is merely a toy example. Or perhaps you had a coursework question to implement least squares regression to a linear model from scratch?

The AudioTools package was introduced in Maple 10 (2005). So it should work in your Maple 12.

restart;
with(AudioTools):

# A[6] 1760.00 Hz
A := Create(duration=2.0, rate=44100,
            x->evalhf(sin(x*2*Pi*1760.00/44100))):
Write(cat(kernelopts(homedir),"/mapleprimes/sampA.wav"), A):

# B[6] 1975.53 Hz
B := Create(duration=2.0, rate=44100,
            x->evalhf(sin(x*2*Pi*1975.53/44100))):
Write(cat(kernelopts(homedir),"/mapleprimes/sampB.wav"), B):

# A[4] 440.00 Hz
A440 := Create(duration=2.0, rate=44100,
            x->evalhf(sin(x*2*Pi*440.00/44100))):
Write(cat(kernelopts(homedir),"/mapleprimes/sampA440.wav"), A440):

As you can see this is very straightforward. Producing, say, B[4] = 493.88 Hz is done similarly.

And of course you could also write a 1-line procedure that took a frequency as an argument and generated the corresponding pure tone.

Using the App Centre link you gave (with its code from the year 2000 using writebytes, convert, and much Library level code) seems unnecessarily much more complicated, and quite possibly slower and less robust.

L1:=[a,b]:

`/`(op(L1));

                 a/b

L2:=[a,b,c,d]:

`/`(op(L2));

                   a
                 -----
                 b c d

It looks like you're trying to compute something simple (like, say, the min or max real root of some cubics). But you have not provided a clear description of your goal. And your worksheet contains weird unreferenced expressions and parameters, as if part of your inputs were copy&pastes.

Trying to deal with symbolic expressions representing roots (and which contain explicit radicals) by calling solve on a cubic with float coefficients will lead to numeric roundoff error upon subsequent floating point evaluation. That includes the case involving substituition for parameter sigma. One consequence is that the purely real root(s) may be hidden behind nonzero imaginary float artefacts. Even the fact that a cubic has at least one real root may be obscured by this unsound technique.

Here's my guess as to what you're trying to accomplish, using fsolve.

In the following attachment I break up the plotting in order to force use of the value of sigma where the discriminant (with respect to a) of the first cubic changes sign. This allows the curves of two of the "roots" to match up.

restart

interface(warnlevel = 0)

R := proc (ee, s, i) local res; options remember, system; if not s::numeric then return ('procname')(args) end if; res := [fsolve(eval(ee, sigma = s))]; if i <= nops(res) and 0 < res[i] then res[i] else Float(undefined) end if end proc

eq := sigma = 1.004210973*a^2+0.1634980569e-2/a

new := (rhs-lhs)(expand(a*eq))

1.004210973*a^3+0.1634980569e-2-a*sigma

cutoff := fsolve(discrim(new, a))

0.2626543768e-1

P1a := plot([R(new, s, 1), R(new, s, 2), R(new, s, 3)], s = -.2 .. cutoff, size = [500, 400], axes = box, color = black); P1b := plot([R(new, s, 1), R(new, s, 2), R(new, s, 3)], s = cutoff+1.0*10^(-Digits+1) .. 1.2, size = [500, 400], axes = box, color = black)

eq2 := sigma = 1.004210973*a^2-0.1634980569e-2/a

new2 := (rhs-lhs)(expand(a*eq2))

1.004210973*a^3-0.1634980569e-2-a*sigma

P2 := plot([R(new2, s, 1), R(new2, s, 2), R(new2, s, 3)], s = -.2 .. 1.2, size = [500, 400], axes = box, color = black)

plots:-display(P1a, P1b, P2)

 

Download cposfsolvediscrim.mw

Here is an earlier attempt using fsolve, which did not force use of the value for sigma where the discriminant changed sign. Notice the hole in the curves where two of the "roots" to not quite match up.
Download cplotposfsolve.mw

Here is an attachment that plots all the roots (not just a>0), using indexed RootOfs.
Download cplot.mw

Here is a related way to obtain the plot of only the positive roots, using the indexed RootOfs. (This is slower than the above way, using fsolve, I believe.) 

Download cplotpos.mw

The very last 2D Input equation in subsection Projektioner is invalid (corrupted, in the sense that the stored base64 encoded .m represents a blob of Typesetting calls that are not validly displayable).

But I believe that the following attachment contains everything else, recovered. Please check.

recovered_ac.mw

(My attachment above was saved in Maple 2019.1. Just in case there's any issue, here also is the original attachment by the OP, with only the offending XML Equation manually removed from the .mw file. raw_edited_ac.mw )

[edited]
Note, sometimes when the GUI has collapsed everything to a single line then formatting can be recovered by doing the main menubar actions View->Sections->Collapse All Sections followed by View->Sections->Expand All Sections. For the OP's corrupted original that doesn't fix everything, as the invalid 2D Input instance makes the GUI go wrong still.

The problematic 2D Input instance is decoded here. errantTS.mw

If I attempt to form a similar invalid 2D Input typeset column Vector (with simply a minus sign in one entry) then the GUI indicates it's invalid but otherwise doesn't run amok. errant2D.mw So I don't know how to reproduce the bug. 

You're second to last line does t:=evalf(t[1]) but before that `t` is not something which can be indexed. So that line makes no sense.

And so what you're passing to fsolve is not a sensible expression.

Your Question is marked as Maple 2019, so I'll use that.

If I load the linkedlist module from a .mla Library Archive (including the one that resides within a workbook), or if I first savelib it to a .mla, then my output from your examples shows the long form linkedlist:-nil provided that I have the typesetting level set to standard.

If I toggle that level back to extended (ie. its original default for new users) then the output shows the short form nil only.

If I simply read the .mpl source then my output shows the short form nil regardless of typesetting level.

Perhaps this is the source of your problem.

nb. You can set the level as a GUI preference -- automatic across sessions. Or you can programmatically set it for an opened worksheet by issuing the command,
   interface('typesetting'=':-extended'):
in its own execution group.

nb. The level extended is the new default for new users. It behaves better than it did back when standard was the default. IMHO it may be generally best to use the default, unless there is some special cause.

I believe that it is not.

First 6 7 8 9 10 11 12 Last Page 8 of 201