Carl Love

Carl Love

25816 Reputation

25 Badges

10 years, 355 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

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

`⇔`(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). That can be done like this

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

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
----------------
R-squared: 0.9941, Adjusted R-squared: 0.9921

[Vector(3, {(1) = 0.5202557353350474e-1, (2) = -1.0044218089795456, (3) = -0.5852877022519349e-1}), 0.47590683831025216e-1]

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
    )
);

[A = 350.156, k__1 = 0.1e-4, Q__mu = HFloat(1659.4444444444443), Q__sigma = HFloat(95.01461875826149), H__mu = HFloat(225.5677777777778), H__sigma = HFloat(8.940258354457344)]

{Q__0 = -2070.65294425443, Q__s = -3562.31910980619, k__f = 0.479613654503028e-4}, {Q__0 = -13158.1287730397, Q__s = 3562.31910980619, k__f = 0.479613654503028e-4}

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
);

HFloat(0.047590683831025646)

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
);

HFloat(33.08811409172967)

NULL

Download SolveBlockQuestion.mw

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))):
addproperty(K, {RealRange(Open(-2), Open(-1))}, {BottomProp});
constants,= K;  #Or constants:= constants, K
protect(K);

 

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.

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.

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]

 

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. 

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.

Suppose that your DataFrame is named DF. To remove all rows for which any cell is undefined, do

DF[andseq(DF[C] <>~ undefined, C= with(DF))];

Great Question, by the way; it's not at all "simplistic". It took me about an hour of experimentation and help-page research to figure out how to do it.

Here's a significantly faster procedure:

OEIS_A219954:= proc(n::posint)
option remember;
local L:= ilog2(n), k:= n-1;
   `if`(n=1, 0, thisproc(k) + `if`(2^L=n, 3^L - n/2, 3^add(Bits:-Split(k))))
end proc
:

Edit: Changed the procedure from an arrow operator to a proc ... end proc to accomodate 2D Input (and I hope that I don't live to regret it) or older Maple.

First, I want to caution you to make functionList a proper Maple list by using square brackets [ ] instead of curly braces { }. The braces make it a set, and you can't control its order. (This is the only reason to make it a list rather than a set.)

The procedure below will take any of your functions and return side-by-side 3D plots of its real and imaginary parts.

ReIm_plot:= proc(f::algebraic, Zr::(name=range(complexcons)))
local z:= lhs(Zr), R:= rhs(Zr), x, y, RI:= [Re,Im];
    plots:-display(
        <
            plot3d~(
                RI(eval(f, z= x+I*y)), (x,y)=~ map~(RI, R)[], 
                labels=~ `[]`~(RI[](z), RI(evaln(:-f)(z))),
                title=~ typeset~(["Real", "Imaginary"], " part of ", f),
                _rest
            )
        >^%T
    )
end proc
:

functionList:= [2 + z, z^2 - 3*z, -z^3 + 4];

You can do one function via

ReIm_plot(functionList[3], z= -1-I..1+I));

Or you can do them all at once via
print~(ReIm_plot~(functionList, z= -1-I..1+I)):  #Note the colon terminator
 

functionList:= [2 + z, z^2 - 3*z, -z^3 + 4]:

ReIm_plot:= proc(f::algebraic, Zr::(name=range(complexcons)))
local z:= lhs(Zr), R:= rhs(Zr), x, y, RI:= [Re,Im];
    plots:-display(
        <
            plot3d~(
                RI(eval(f, z= x+I*y)), (x,y)=~ map~(RI, R)[],
                labels=~ `[]`~(RI[](z), RI(evaln(:-f)(z))),
                title=~ typeset~(["Real", "Imaginary"], " part of ", f),
                _rest
            )
        >^%T
    )
end proc
:

print~(ReIm_plot~(functionList, z= -1-I..1+I)):

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Download ReIm_plot.mw

I think that you want cylindrical coordinates. But showing the discontinuity (division by 0) elegantly is more difficult in 3d than 2d, regardless of your coordinate system.

plot3d(
    [r, theta, 1/(r^2*sin(theta)^2)], r= 0..2, theta= -Pi..Pi,
    coords= cylindrical, shading= zhue, view= [-2..2, -2..2, 0..99]
);