mmcdara

4029 Reputation

17 Badges

6 years, 183 days

MaplePrimes Activity


These are answers submitted by mmcdara



The position of the legend and the line spacings can be adjusted as you want.
Hint : 

  • op~(2, [plottools:-getdata(p1)]) gets the x and y ranges of all the plots
  • op~(2, op~(2, [plottools:-getdata(p1)])) gets the y ranges of all the plots
     

restart:

with(plots):

equ1 := BesselJ(sqrt(17)/2, 10*sqrt(t)*sqrt(2))/t^(1/4) + BesselY(sqrt(17)/2, 10*sqrt(t)*sqrt(2))/t^(1/4):

equ2 := BesselJ(sqrt(17)/2, 10*sqrt(t)*sqrt(2))/t^(1/4) + 5*BesselY(sqrt(17)/2, 10*sqrt(t)*sqrt(2))/t^(1/4):

equ3 := BesselJ(sqrt(17)/2, 10*sqrt(t))/t^(1/4) + 5*BesselY(sqrt(17)/2, 10*sqrt(t))/t^(1/4):

tmax   := 30:
colors := ["Red", "Violet", "Blue"]:

p1 := plot([equ1, equ2, equ3], t = 0 .. tmax, labels = [t, T[2](t)], tickmarks = [0, 0], labelfont = [TIMES, ITALIC, 12], axes = boxed, color = colors):

ymin := min(op~(1, op~(2, op~(2, [plottools:-getdata(p1)])))):
ymax := max(op~(2, op~(2, op~(2, [plottools:-getdata(p1)])))):
dy   := 2*ymax:

legend1 := typeset(C[3] = 1, ` , `, C[4] = 1, ` , `, Omega^2 = 50):
legend2 := typeset(C[3] = 1, ` , `, C[4] = 5, ` , `, Omega^2 = 50):
legend3 := typeset(C[3] = 1, ` , `, C[4] = 5, ` , `, Omega^2 = 25):

p2 := seq(textplot([tmax-2, ymax-k*dy/20, legend||k], align=left), k=1..3):

p3 := seq(plot([[tmax-2, ymax-k*dy/20], [tmax-1, ymax-k*dy/20]], color=colors[k]), k=1..3):
display(p1, p2, p3, view=[default, -ymax..ymax], size=[800, 500])

 

 


 

Download Legend_Inside.mw



Leap year not accounted for (every february here has 28 days)
 

restart:
month := [seq(cat("0", k), k=1..9), seq(cat("1", k), k=0..2)]:
day   := convert~([$1..9, seq(seq(cat(i, k), k=1..9), i=1..2), 30, 31], string):
bounds := [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]:
alias(CR=combinat:-randperm):
m := CR(month)[1];
d := CR(day[1..parse(m)])[1];
y := parse(StringTools:-Reverse(cat(d, m)));
                              "09"
                              "7"
                              907
cat(cat(d, "/", m), "/", y)
                           "7/09/907"

 


For non linear fit I usually prefer using Optimization:-NLPSolve instead of Statistics:-NonlinearFit.
As a rule Optimization:-NLPSolve is much faster than Statistics:-NonlinearFit and, as it is very flexible to define constraints, it avoids doing this loop on n.

Using  Optimization:-NLPSolve will enable you to select easily the "correct best solutions" (in the second attempt I used c >=-100).

I let you compare the execution time using  Optimization:-NLPSolve or your code (in black at the end of the file).

 

restart;

with(Statistics):

j := Matrix(20, 3, [[70, 8, 8.468140006], [70, 10, 4.105515432], [70, 12, 2.36261199], [70, 14, 1.422093923], [70, 16, 0.9], [100, 8, 20.47249229], [100, 10, 9.618450629], [100, 12, 5.360869165], [100, 14, 3.399312905], [100, 16, 2.640399788], [130, 8, 35.90466304], [130, 10, 17.62958097], [130, 12, 9.828362586], [130, 14, 5.866863694], [130, 16, 3.799262645], [160, 8, 57.31814648], [160, 10, 34.49692774], [160, 12, 15.39340528], [160, 14, 9.991012951], [160, 16, 6.049343013]]):

# last command killed after 30 seconds
# eq := NonlinearFit(u*h^n*exp(-b*d)/(-c*d + h), j, [h, d], output = [parametervalues, residualsumofsquares]);

# define the model f(h, d) and the objective function obj (sum of squared residuals)

f   := (h, d) -> u*h^n*exp(-b*d)/(-c*d + h):
obj := add( (f(j[k, 1], j[k, 2]) - j[k, 3])^2, k=1..numelems(j[..,1]) ):

# Use Optimization:-NLPSolve instead of Statistics:-NonlnearFit (the former is more versatile)

opt := Optimization:-NLPSolve(obj, {n>=2, n<=3}, method=sqp, iterationlimit=1000)

[21.6509289224931862, [b = HFloat(0.2028991836842215), c = HFloat(-6477.6480267378065), n = HFloat(2.3544626835278906), u = HFloat(99.32223640618867)]]

(1)

BestModel := eval(f(h, d), opt[2])

HFloat(99.32223640618867)*h^HFloat(2.3544626835278906)*exp(-HFloat(0.2028991836842215)*d)/(HFloat(6477.6480267378065)*d+h)

(2)

plots:-display(
  ScatterPlot3D(j, symbol=solidsphere, color=green, symbolsize=15),
  plot3d(BestModel, h=(min..max)(j[..,1]), d=(min..max)(j[..,2]), style=wireframe, color=blue)
);

 

# add an extra constraint on c

opt := Optimization:-NLPSolve(obj, {n>=2, n<=3, c >=-100}, method=sqp, iterationlimit=1000);
BestModel := eval(f(h, d), opt[2]):
plots:-display(
  ScatterPlot3D(j, symbol=solidsphere, color=green, symbolsize=15),
  plot3d(BestModel, h=(min..max)(j[..,1]), d=(min..max)(j[..,2]), style=wireframe, color=blue)
);

[20.2914105165715100, [b = HFloat(0.24364537500713906), c = HFloat(-20.866765615558734), n = HFloat(2.7593674638195393), u = HFloat(0.11018019031068993)]]

 

 

restart;

with(Statistics):

j := Matrix(20, 3, [[70, 8, 8.468140006], [70, 10, 4.105515432], [70, 12, 2.36261199], [70, 14, 1.422093923], [70, 16, 0.9], [100, 8, 20.47249229], [100, 10, 9.618450629], [100, 12, 5.360869165], [100, 14, 3.399312905], [100, 16, 2.640399788], [130, 8, 35.90466304], [130, 10, 17.62958097], [130, 12, 9.828362586], [130, 14, 5.866863694], [130, 16, 3.799262645], [160, 8, 57.31814648], [160, 10, 34.49692774], [160, 12, 15.39340528], [160, 14, 9.991012951], [160, 16, 6.049343013]]):

minRSS := +infinity:
for n from 2 by 0.001 to 3 do

    eq := NonlinearFit(u*h^n*exp(-b*d)/(-c*d + h), j, [h, d], initialvalues = [b = 0.1, c = -30, u = 0.1], output = [parametervalues, residualsumofsquares]);

    if eq[2] < minRSS then

        #print(n, eq[2], eq[1]);
        optRSS := [''n''=n, ''minRSS'' = eq[2], eq[1]]:
        minRSS := eq[2]:

    end if;

end do:

optRSS;

[n = 2.759, minRSS = 20.29141110, [b = HFloat(0.24360946347455648), c = HFloat(-20.898174204286263), u = HFloat(0.11043917832727822)]]

(3)

 


 

Download ConstrainedNonLinearFit.mw

 

To complete @tomleslie 's answer

If your matrix is that big, save it in an m file (binary file): you will spare time and bytes.

Let M the name of the matrix and f := "//../MyFile.m" some m file

  • save M, f   will save the matrix

To read it in a new worksheet

  • define f := "//../MyFile.m"
  • read f  will read the content of f
    This command automatically creates a matrix whose name is the same than the name used in the save command.
     

To verify this :

restart:
f := "//..../MyFyle.m";
read f;          # no echo of what is done
anames(user);    # will display the user variables, here f and M
                 # Useful if you don't remember the name the matrix had in the initial worksheet
                 # or if you share your work with other people

 


Not perfect but I have only Maple 2015 right now.
(advice: browse the questions to seek the works @acer has already done on related problems)

Download Colors.mw


It still remains things that only you can do correctly (for instance what D(T)^2 is supposed to mean???)

Here is a file with syntax errors fixed.

restart;

A := 7.17;

B := 2.56*10^(-3);

C := 0.08*10^5;

_local(D): # not a MAPLE's syntax
local D:

D := 0*10^(-6);

T_0 := 298;

G0 := -71.398;

S0 := 45.106;

Hf := -57.95;

cp := A + B*T + C/T^2 + D(T)^2;   # WRONG
cp := A + B*T + C/T^2 + D*T^2;    # Is that what you mean?
                                  # If not write yourself the term involving D

7.17

 

0.2560000000e-2

 

8000.00

 

0

 

298

 

-71.398

 

45.106

 

-57.95

 

7.17+0.2560000000e-2*T+8000.00/T^2

 

7.17+0.2560000000e-2*T+8000.00/T^2

(1)

# the integration variable can't be an integration bound

`&Delta;H` := T -> int(cp, z = T_0 .. T):
`&Delta;S` := T -> int(cp/z, z = T_0 .. T):

G := T -> -S0*T - T*`&Delta;S`(T) + Hf + `&Delta;H`(T):


# in order to plot something TMAX must be given a numerical value first.

TMAX := 400:
p1 := plot(G(T), T = T_0 .. TMAX, color=red, legend=typeset(''int(cp/z, z = T_0 .. T)'')):

# But maybe you wanted to write `&Delta;S` this way?

`&Delta;S` := T -> int(cp/T, z = T_0 .. T):

p2 := plot(G(T), T = T_0 .. TMAX, color=blue, legend=typeset(''int(cp/T, z = T_0 .. T)'')):

plots:-display(p1, p2, gridlines=true)

 

 

Download Errors_fixed.mw

Are you familiar with linear regression in Maple?
If not here is a simple example as a starting point.
LinearFit_Example.mw

You will see in the first one that the requirement R^2 > 0.9 is not met (while it is in the second one).
A way to improve its value is to complexify the model by adding quadratic regressors (x[i]^2) or 2nd order interactions among regressors (x[i]*x[j], j<>i), or both of them.

As you have 32 variables there exist (32*31)/2 = 496 terms of the form  x[i]^2  and x[i]*x[j] (j<>i).
If the linear model without constant term depends on 32 parameters, the complete quadratic model depends on 528 parameters.
Note that this implies the number of observations is at least equal to 528 for the Least Squares method to work.

Remark: this constraint can be removed for come specific regression methods (see for instance  Partial_least_squares_regression)

The number of parameters a linear model depends upon is often named the "complexity of the model".
Some advanced regression methods (Stepwise_regressionRidge_regressionLasso_(statistics) and a little more) are aimed to seek for low complexity models which have nevertheless a quality of approximation (R^2 for instance) close to a far more comlex model (for instance, in your case a model of complexity 63 could be almost as good as the complete quadratic model of complexity 528).

Nevertheless, be aware that the closest to the number N of observations the number of parameters P is, the higher the value of R^2.
For P=N one has even R^2 = 1!!!

If you want to go further let me know. 

This comple tomleslie's reply

51 linear equations is not that big for Maple....
... but the output will be very large (particularly if you use w's expression and if you use integers or rationals instead of floats
So the first thing to do is not to define w! (this will be done after simplifying the solution in an ad hoc way).

Look to this attached file (no particular precaution taken to enhance the performances)

Download example.mw

  1. Search Assistant in the help pages
  2. Select the item worksheet, documenting, assistantsmenu (assistant)
  3. Select Back-Solver
  4. In the Back-Solving Assistant worksheet
    1. copy one the three (I guess) formulas that are given
    2. paste it in the empty field above the button Proceed to Back-Solver
  5. Click on the button Proceed to Back-Solver, this will open a new worksheet named 'untitled..."
  6. Jump into this latter, Enjoy

A suggestion: before running your computation, open the task monitor and be ready to kill the application (or the mserver process) before  it all goes wrong (in some cases using the "stop" button doesn't work).

Computing the determinant of a 3x3 matrix should be an easy task, unless each element of the matrix has a complex and lengthy expression.You seem to say that is your case? If it is so, unless some simplifications occur, the determinant will be even more complex. It could even take pages and pages to be displayed?
In such a case, what is the point to compute its formal expression if you are not able to understand it?

So another suggestion: Do the operations step by step, for instance simplify a22*a33, simplify a23*a32, form their difference, simplify (or collect/factor/normal(ize)), multiply by a11 ; do the same for the other minors
For instance

restart:

# Let p a list/vector of arbitrary length.
# Let's suppose your matrix is something like that

M := Matrix(3$2, (i, j) -> f[i,j](p));

M := Matrix(3, 3, {(1, 1) = f[1, 1](p), (1, 2) = f[1, 2](p), (1, 3) = f[1, 3](p), (2, 1) = f[2, 1](p), (2, 2) = f[2, 2](p), (2, 3) = f[2, 3](p), (3, 1) = f[3, 1](p), (3, 2) = f[3, 2](p), (3, 3) = f[3, 3](p)})

(1)

# Let A a "simple" symbolic matrix and detA is determinant

A    := Matrix(3$2, symbol=a):
detA := LinearAlgebra:-Determinant(A);

a[1, 1]*a[2, 2]*a[3, 3]-a[1, 1]*a[2, 3]*a[3, 2]-a[1, 2]*a[2, 1]*a[3, 3]+a[1, 2]*a[2, 3]*a[3, 1]+a[1, 3]*a[2, 1]*a[3, 2]-a[1, 3]*a[2, 2]*a[3, 1]

(2)

# rewritting rules

RR := [entries(A=~M, nolist)]

[a[1, 1] = f[1, 1](p), a[2, 1] = f[2, 1](p), a[3, 1] = f[3, 1](p), a[1, 2] = f[1, 2](p), a[2, 2] = f[2, 2](p), a[3, 2] = f[3, 2](p), a[1, 3] = f[1, 3](p), a[2, 3] = f[2, 3](p), a[3, 3] = f[3, 3](p)]

(3)

# obviously the determinant detM of M is

detM := 'eval(detA, RR)'

eval(detA, RR)

(4)

# Now let's suppose you want to study the behaviour of detM when all
# the values od p but the nth are given a fixed value "val".
# To do this evaluate RR for p[i]=val[i], i:1..n-1, n+1.. and
# compute eval(detA, eval(RR, ...))

# Meanwhile you can evaluate separately the 6 terms of detA and try
# to simplify them one way or another.
# For instance (use simplify/collect/factor, ...)

simplify(eval(op(1, detA), RR));

# and keep doing this for all the terms.
# If you are lucky enough you could obtain something tractable and understandable

 

 

Download determinant.mw

Create an alias?

 

restart:
psi;
                              psi
alias(psi=psi(t))
                              psi

# psi can be us without explicit reference to t:
diff(psi, t)
                             d     
                            --- psi
                             dt    
dsolve(diff(psi, t)=t)
                              1  2      
                        psi = - t  + _C1
                              2         

 

To project a function  f(x,y) defined over some (cartesian) domain [a, b]x[c, d] you need before all to have some mapping from this plane domain onto the sphere S2.
As this is mathematically impossible you can only layer a structure on S2 in an approximate way.
 
A simple case is  [a, b]x[c, d] = [0, 2*Pi]x[-Pi, Pi] (each bounded rectangular plane domain can be mapped onto this one).
Generally this "layering" introduces singularities. To avoid them f has to be x-periodic, y-periodic, and both f(x, -Pi) and f(x, Pi) must be constant (unique values at the north and south poles);

Here is a way to "project"  a function such that f(x,y) on the unit sphere S2

projection_onto_S2.mw



 

that I found elegant but which has the drawback of becoming dramatically slow as soon as N is "slightly large".
So this is just for fun

restart:
with(PolyhedralSets):
with(ExampleSets):
N := 3:
H := Hypercube(N):
V := convert(VerticesAndRays(H)[1], Matrix)

In the command

plots[implicitplot](something, P=1..4, Q=0..3, ...);

something is aimed to be an expression which depends continuously on P and Qn (or at least which is continuous wrt P and Q but on set of null measure).
In your case this something is almost everywhere equal to 

1/(1+1.5*exp(0.01*s))

When P varies from 1 to 4 and Q varies from 0 to 3 the couple (P, Q) will take a lot of values whose only (1, 0), (1, 3), (4, 0) and (4, 3) are integers.
Thus  the definition FC:= proc(s, P::integer, Q::integer) is almost surely going to generate the error message you got.
One way coud br to write FC:= proc(s, P, Q) instead.
But in this case what are you expecting to plot?

The command 

plots[implicitplot](something, P=1..4, Q=0..3, ...);

will be equivalent to 

plots[implicitplot](1/(1+1.5*exp(0.01*s[i])), t=, P=1..4, Q=0..3, ...);

for some given value of s.
Which means you are trying to implicitplot a constant wrt P and Q.
Look for instance what happens wifh this

plots:-implicitplot(1,  P=1..4, Q=0..3)


I don not understand what you are trying do plot, but I wonder if this isn't what you would want... by any chance

restart

M := table([
        And(P=2,Q=0) = 1/(1 + 2.1*exp(0.01*s)),
        And(P=2,Q=1) = 1/(2 + 1.5*exp(0.0105*s)),
        And(P=3,Q=0) = 1/(2 + 3.2*exp(0.011*s)),
        And(P=3,Q=1) = 1/(2.5 + 2*exp(0.012*s)),
        And(P=3,Q=2) = 1/(3 + 1.2*exp(0.011*s)),
        And(P=4,Q=0) = 1/(3 + 10*exp(0.012*s)),
        And(P=4,Q=1) = 1/(2 + 5*exp(0.011*s)),
        And(P=4,Q=2) = 1/(1.5 + 3*exp(0.011*s)),
        And(P=4,Q=3) = -33000/(-150000 - 60000*exp(0.011*s)),
        others       = 1/(1+1.5*exp(0.01*s))
     ]):  

conds := sort([indices(M, nolist)])

[others, And(P = 2, Q = 0), And(P = 2, Q = 1), And(P = 3, Q = 0), And(P = 3, Q = 1), And(P = 3, Q = 2), And(P = 4, Q = 0), And(P = 4, Q = 1), And(P = 4, Q = 2), And(P = 4, Q = 3)]

(1)

seq(M[c], c in conds)

1/(1+1.5*exp(0.1e-1*s)), 1/(1+2.1*exp(0.1e-1*s)), 1/(2+1.5*exp(0.105e-1*s)), 1/(2+3.2*exp(0.11e-1*s)), 1/(2.5+2*exp(0.12e-1*s)), 1/(3+1.2*exp(0.11e-1*s)), 1/(3+10*exp(0.12e-1*s)), 1/(2+5*exp(0.11e-1*s)), 1/(1.5+3*exp(0.11e-1*s)), -33000/(-150000-60000*exp(0.11e-1*s))

(2)

r := rand(0. .. 1.):

plots:-display(
  seq(plot(M[c], s=0..100, color=ColorTools:-Color([r(), r(), r()]), legend=c), c in conds),
  legendstyle=[location=right],
  size=[500,400]
)

 

 

NULL

Download proposal.mw

In case you would like to order the rows in such a way that only one element changes from a row to the next one you can use this 

restart
f := proc(N)
  local b := convert~(combinat:-graycode(N), binary):
  local L := [ seq(parse~(StringTools:-Explode(sprintf(cat(" %.3d"), b[k]))), k=1..2^N) ]:
  local M := convert(L, Matrix);
  return 2*~M -~ 1
end proc:

f(1), f(2), f(3)

 

3 4 5 6 7 8 9 Last Page 5 of 33