mmcdara

7891 Reputation

22 Badges

9 years, 55 days

MaplePrimes Activity


These are answers submitted by mmcdara

As tomleslie opened the gate, here is a model which is infered from the observation of the data and which is aimed to respect some of their feautures:

  • Y is (seems to be) positive
  • Y is monotonically increasing (unless if Y[298]=Y[297]-1-

More of this a close examination of your data leads me to infer there are two adjacent X ranges:

  • for X fro 1 to (about) 220, Y evolves in some way
  • for X > 220 the evolution of Y seems to be the same as in the previous range, maybe just shifted and scaled in a propper way

A model that respect these requirements and observations is

th  := x -> a__1*(tanh((x-c__1)*b__1)+1) + a__2*(tanh((x-c__2)*b__2)+1):

If you discard the positivity constraint you can use this more general model

th := x -> a__1*(tanh((x-c__1)*b__1)+d__1) + a__2*(tanh((x-c__2)*b__2)+d__2)

Here is the result for this latter


Fit_of_your_data.mw

Other solutions (better or worse) can be obtained by using similar models:

  • replace tanh by arctan
  • replace tanh by a sigmoid function

As a rule "before using a nonlinear fir, always look if the model is not linear under a proper transformation".
This will potentialy avoid numerical problems.
In your case the transformation is obvious:

y = A*exp(r*x) --> log(y) = log(A) + r*x

Thus:
 

beta := LinearFit([1, x], X, log~(Y), x, output=parametervalues);

Model := x -> exp(beta.<1, x>):
expand(Model(x));

    117.5619317177743 * exp(0.02875273320702703 * x)

plots:-display(
  ScatterPlot(X, Y), 
  plot(Model(x), x=(min..max)(X), legend=evalf[4](expand(Model(x)))), 
  view=[default, 0..Model(X[-1])]
);




The fit is nevertheless rather poor.
So let's use this fit as a starting point and look if one can do better?
IMO Fit is not very efficient and I prefer doing the things "by hand" and to use
Optimization:-Minimize instead:

Prior := x -> A*exp(x*r):

# Objective function 

J := add((Y - Prior~(X))^~2):

opt := Optimization:-Minimize(J, initialpoint = {A=110, r=0.03});
                                                 
  [3.26776e10  , [A = 10049.7774, r = 0.00735]]

# For information, this is the value of J for Model(x)
J := add((Y - Model~(X))^~2);
                  HFloat(4.49e12)

Post := eval(Prior(x), opt[2]);

plots:-display(
  ScatterPlot(X, Y), 
  plot(Post, x=(min..max)(X), legend=evalf[4](expand(Model(x)))), 
  view=[default, 0..Model(X[-1])]
)



LinearFit.mw

Even if Post is better than Model for it gives an objective function 100 times smaller, I'm not sure that we can't do even better by choosing different initial values for A and r...


For an unknown reason DrawGraph transforms each weight w as convert(evalf(w), string) and generates a TEXT structure for each weight which looks like this one

# Here is how the weight on the edge from not A to C is displayed
# Note that 1/7 has been transformed as "0.143"
op(-5, P)
    TEXT([0.7752048164, 0.2971539501], "0.143", ALIGNRIGHT, ALIGNBELOW, FONT(HELVETICA, 11))

whereas what is expected is

TEXT([0.7752048164, 0.2971539501], 1/7, ALIGNRIGHT, ALIGNBELOW, FONT(HELVETICA, 11))


Knowing this is quite easy to recover the initial form

# construct the graph WITHOUT weights displayed 

P__unweighted  := [ op(DrawGraph(G, style=tree, root=Omega, showweights=false)) ]:

# Extract the part of the plot that concerns the weights

EdgeWeights    := map(u -> if `not`(u in P__unweighted) then u end if, [ op(P) ]):

# Recode the weights as rationals

NewEdgeWeights := map( u -> subsop(2=convert(parse(op(2, u)), rational, 3), u), EdgeWeights):

# Reassemble the "weighted" structure

P__weighted    := P__unweighted[], NewEdgeWeights[]:

# Plot it
# as in my answer to your initial question I syill use Maple 2015, 
# so the rendering might be slightly different from yours

HighlightVertex(G, Vertices(G), white);
plots:-display(R(PLOT(P__weighted)), size=[600, 600]);

 
PS: I suspect the rational to float transformations are justified to avoid complicating the algorithm wich calculates vertex positions.
(the weight information "takes" on line instead of 2).
In my opinion a transformaton of the form 

r := 1/7:
rr := numer(r) %/ denom(r)

whould be better

PS: As I said you while answering your first question about probability graphs, you cann replace a vertex name like A by defining  A:=`#mo(A)` in order to print all the vertices with the same font. 

Probability_Graph_2.mw

I understand from your previous question that you are familiar with LaTeX?
In this case you could be find a few ideas here

LateX_output.mw

Basically I generate the LaTeX code of a "tabular" structure directly in Maple.
The only prerequisite is yo know LaTeX.
I'm responsible of a computational code written in Maple and I have several test cases I need to run periodically to insure there is no code regression. I use to use the trick described in the attached file to generate automatically "execution reports" and insure traceability.

PS: I did not compile the tex cource to verify I didn' make mistakes, but I think I'm right

Hi, 

It's a consequence of the (rather comple) algorithm implicitplot uses (type showstat(`plots/implicitplot`)  in a worksheet and look at it).
In short, if numpoints if specified, implicitplot begins by computing a grid whose dimension will be nx by ny, where 

nx, ny := `$`(isqrt(numpoints+3)+1,2)

(for numpoints=1000 this gives a 33 by 33 grid).
This is the grid the curve is interpolated on.
The resolution intervenes later to determine the number of point where F(x, y) is to be evaluated.

For gridrefine=3 F(x,y) is calculated at 515 points, while 4119 points are used for gridrefine=6.
(Note these numbers are twice the numbers I obtained by examining the output of showstat(`plots/implicitplot`) ).

Now there is a third option named resolution, wich plays a role when explitely set to a value its default value is 0).
By setting resolution = k times 'the number of points where F(x,y) is computed' (k integer > 0) you can force implicitplot to really use all the points (and more) to get a ... better resolution of the curve.

In the attached file:

  • p3 is the plot obtained with gridrefine=3,
  • p6_0 is the plot obtained with gridrefine=6, and resolution=0 (default value),
  • p6_4 is the plot obtained with gridrefine=6, and resolution=4*4119 (or about).

    nonsmooth_implicitplot.mw

A remark : 

M[i, j] := M[i, j] + something

# is useless, it's enough to write

M[i, j] := something


This is a shorter way to obtain the matrix M:
(maybe "expr" is less readable)
 

restart;


A := Matrix(2, 3, [[4, x, 1/4*x^2], [2, 3, x]]);
B := Matrix(3, 2, [[3, 5], [x, 2], [1/2*x^2, -x]]);

a := LinearAlgebra[ColumnDimension](A):
b := LinearAlgebra[ColumnDimension](B):

A := unapply(A, x):
B := unapply(B, x):

alpha := k -> x -> A(x)[k, ..]:
beta  := k -> x -> B(x)[.., k]:


expr := (i,j) -> alpha(1)(1/2)[i] * beta(1)(1/2)[j]
                 +
                 alpha(2)(1/2)[i] * beta(2)(1/2)[j]
                 +
                 int(alpha(1)(x+1/2)[i] * beta(2)(1/2*x)[j], x =-1/2..1/2);

M := Matrix(a$2, (i, j) -> expr(j, i));

AnotherWay.mw

The first example can be proved his way

with(Student[Precalculus]):
CompleteSquare(a^2+2*a*b+b^2, a);
CompleteSquare(a^2+2*a*b+b^2, b);
                                   2
                            (b + a) 
                                   2
                            (b + a) 


Be patient for the second...

It works correctly, maybe you have made a false manipulation

restart:
interface(version)
Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895

diff(y(x), x, x) = 1/y(x) - x*diff(y(x), x)/y(x)^2;
sol := dsolve({%, y(0)=1, D(y)(0)=0}, numeric);
plots:-odeplot(sol, [x, y(x)], x=0..1)


ode.mw


By the way, I've given a reply yesterday to your old unanswerd problem Trying to make a simulated orbit accurate for thousands of orbits, have you seen it?

Concerning your current problem:

  • error 1: "omega(t)" is not defined before the plot (could it be the solution of your first problem?)
  • error 2 : the lines
    thetas := eval(`mod`(theta(t), 2*Pi), soln); 
    omegas := eval(omega(t)/omegaH, soln)
    should be replaced, for instance by 
    omegas := s -> eval(omega(t)/omegaH, soln(s));
    pp := 2*Pi:
    thetas  :=  s -> eval(theta(t)-pp*round(theta(t)/pp), soln(s));
    Further, lines 
    thetavalues[n] := fsolve(thetas(t) = n*TH, t); 
    omegavalues[n] := fsolve(omegas(t) = n*TH, t);

    must be changed by

  • thetavalues[n] := fsolve('thetas(t)' = n*TH, t);
    omegavalues[n] := fsolve('omegas(t)' = n*TH, t)
    At last, I think you should use the  code contained in my reply to Trying to make a simulated orbit accurate for thousands of orbits ; it seems you are doing the same things and I've fixed the problem you then faced
  • errors in "Investigation 1": what did you want to do with these characters (* and *) 
  •  

I put the dsolve and display commands within an "if... end if" structure to avoid generating a huge file that cannot even be loaded on this site. 
There are minor modifications compariring to your original file

  • I have set Digits to 15 at the top of the worksheet
    (I'm not sure this really help improving the results)
  • I have defined T=100*TH BEFORE the call to dsolve
  • As T is now defined, I cab use the option range=0..T within dsolve.
    (at the same time I have replaced rkf45 by rosenbrock, but the improvement, if any, is extremely small)
  • In the display commands:
    • I've used the option refine=1 (be carefull, refine=1 generates a matrix of size 572000 by 2 and refine=2 generates anerror probably do to an even larger matrix).
    • I've splitted the range 0..T into 10 consecutives ranges of equal length T/10

Project_reply.mw

 

Hi, 

There exists no build-in procedure to do this, but you can use this code (done in Maple 2015, a lot of improvements for customizing a graph exist in earlier versions, so I'm not going to spend time for a better rendering)

 

restart
interface(version)
Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 

   21 2015 Build ID 1097895

with(GraphTheory):

R := plottools:-transform((x,y) -> [-y,x]):

A    := `#mo(A)`:  # to have the same font aas notA
B    := `#mo(B)`:
C    := `#mo(C)`:
notA := `#mover(mo(A),mo("&#95;"))`:
notB := `#mover(mo(B),mo("&#95;"))`:
notC := `#mover(mo(C),mo("&#95;"))`:

# Step 1
#
# A probability graph is either represented as:
#   1.  a weighted directed graph (Markov processes point of view)
#   2.  a weighted undirected graph (Probability point of view, but in this 
#       case there are implicit rules to avoid representing edges as oriented
#       arcs
#
# Set the boolean variable I_Prefer_Directed_Graphs as you want depending
# on your preference.

I_Prefer_Directed_Graphs := true:  # false

if I_Prefer_Directed_Graphs then

  # Step 1

  DG := Graph({
                [[Omega, A], 1/4], [[Omega, notA], 3/4], 
                [[A, B], 1/2], [[A, notB], 1/2], 
                [[notA, C], 0.237], [[notA, notC], 1-0.237]
             }):

  # Step2 
  #
  # But only undirected graphs can be represented in tree style,
  # So I create here the 

  G := Graph(convert~(Edges(DG), set)):
  P := DrawGraph(G, style=tree, root=Omega):

  # Step 3
  #
  # Now I get the vertices position of the tree graph G and set them as
  # vertices positions of the graph DG

  V  := GetVertexPositions(G):
  GD := Graph(convert~(Edges(G), list)):
  SetVertexPositions(DG, V):
  HighlightVertex(DG, Vertices(DG), white);

  DP := DrawGraph(DG):
  print( plots:-display(R(DP), size=[600, 600]) );

else
  G := Graph({
               [{Omega, A}, 1/4], [{Omega, notA}, 3/4], 
               [{A, B}, 1/2], [{A, notB}, 1/2], 
               [{notA, C}, 0.237], [{notA, notC}, 1-0.237]
            }):
  HighlightVertex(G, Vertices(G), white);
  P := DrawGraph(G, style=tree, root=Omega):
  HighlightVertex(G, Vertices(G), white);
  print( plots:-display(R(P), size=[600, 600]) );

end if:


Probability_Graph.mw

You might find useful informations here
How to redirect showstat(....) into a file ?

does this suit you?

restart:
expr:=y^2+y^3*sin(x)+3*x*y^5;
convert(expr, horner, y);
#
expr:=y^2*f(y)+y^3*sin(x)*g(y)+3*x;
convert(expr, horner, y);

horner.mw

Maybe this?
 

restart
assume(x::real, y::real):
z := x+I*y:
f := abs(z-1)+abs(z-I)-2:
plots:-implicitplot(f, x = -2..2, y= -2..2): # OK
g := abs(z-I)-abs(z+1)-2:
plots:-implicitplot(g, x = -2..2, y= -2..2): # KO
sol := [solve(g, y)];
plots:-complexplot(sol, x=-10..10, color=[red, blue]): # OK

Download Z.mw

2*sin(x)*cos(x) = sin(2*x)
next 
sin(2*x) = 2*sin(x)*cos(x)

# 1/  2*sin(x)*cos(x) = sin(2*x)

combine(2*sin(x)*cos(x), trig)

sin(2*x)

(1)

# 2/  sin(2*x) = 2*sin(x)*cos(x)
#     write sin(2*x) as sin(x+y), next set y=x

expand(sin(x+y));
eval(%, y=x)

sin(x)*cos(y)+cos(x)*sin(y)

 

2*sin(x)*cos(x)

(2)

 


 

Download prove.mw

First 47 48 49 50 51 52 53 Last Page 49 of 65