mmcdara

6149 Reputation

17 Badges

9 years, 62 days

MaplePrimes Activity


These are answers submitted by mmcdara

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("_"))`:
notB := `#mover(mo(B),mo("_"))`:
notC := `#mover(mo(C),mo("_"))`:

# 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

Some typos in your pseudo code... here corrected
 

with(plots):

with(plottools):

f:=(x-1)/(x+2):

p1 := plot(f,x=-3..3):

p2 := circle([0,0],1,color=blue):

display(p1,p2, scaling=constrained, view=[default, -5..5]);

 

 


 

Download SomeTypos.mw

A possibility is to load the file manually and then read it (look to packages FileTools or ExcelTools, or use readdata).
a smarter one is to use the package URL :

URL:-Get( "https://www.gw-openscience.org/GW150914data/P150914/fig1-observed-H.txt" ):

... but this command returns an uncomprehensible error (maybe it will wotk for you?)
Error, (in URL:-Get) unknown SSL protocol error in connection to www.gw-openscience.org:-9836
 

There is no error with Maple 2015 (and I think with more recent versions neither)
You will find in the attached file 3 different syntaxes:

  • with curly brackets,
  • with square brackets
  • without any brackets.

All of them work correctly but the result may depend on the syntax and the order ( (k, k1) vs (k1, k) ) you use


Maple_2015.mw



 

A way:


 

restart

interface(version)

`Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895`

(1)

I__D__moy := -(-1 + R)*I__L__moy

-(-1+R)*I__L__moy

(2)

collect(expand(I__D__moy), I__L__moy):
map(sort, %, R, ascending)

(1-R)*I__L__moy

(3)

 


 

Download arrange.mw

Hi,

There are a lot of errors and I doubt you even tried to run this in a Maple worksheet?
Anyway, here is a minimum code to help you do something by yourself:

  • change the values of the parameters as you want
  • define your own plot (odeplot3 is not a procedure in Maple) and customize it as you want

The rest is on you

restart

dsys:={diff(n(t),t)=-2*(n(t)-N)+(u(t)-abs(M)*exp(I*phi))*exp(-2*I*omegap*t/lambda)+conjugate((u(t)-abs(M)*exp(I*phi))*exp(-2*I*omegap*t/lambda)),diff(u(t),t)=-2*(1-I*delta)*u(t)+2*(n(t)-N)*exp(-2*I*omegap*t/lambda)+2*abs(M)*exp(I*phi)}:

res1:=dsolve(dsys union {n(0)=0,u(0)=0},numeric,output=listprocedure);

Warning, The use of global variables in numerical ODE problems is deprecated, and will be removed in a future release. Use the 'parameters' argument instead (see ?dsolve,numeric,parameters)

 

# Maple con solve this system numericaly for some quantities are not defined:

undefined_quantities := indets(dsys, name)

{M, N, delta, lambda, omegap, phi, t}

(1)

# first thing to do is to extract the parameters

params := convert(undefined_quantities minus {t}, list)

[M, N, delta, lambda, omegap, phi]

(2)

res1 := dsolve(dsys union {n(0)=0,u(0)=0}, numeric, parameters=params);

# example of instanciated parameters

data := params=~[1, 1, 1, 1, 1, 1]:
res1(parameters=data):
res1(1)

[t = 1., n(t) = .690491719445859+0.*I, u(t) = -.363285887100282+.879463452206834*I]

(3)

# I don(t know what you want to plot, all the more that n(t) and u(t) are potentially complex

plots:-odeplot(res1, [t, Re(n(t)), Re(u(t))], t=0..1)

 

 


 

Download Try_This.mw


I believe one can print smart tables with Maple too:

(change the format of the column titles to see what happens)

This has been obtained with Maple 2015 and I think it will not work with Maple 2019 or later because the last command 

DocumentTools:-InsertContent(Worksheet(Group(Input( W )))):

has slightly changed from 2015 in more recent versions (unfortunately I do not remember what this change is)
Someone alse here (@acer for instancep will be able to fix this if there is indeed an issue.

 

restart:

with(DocumentTools:-Layout):
with(ScientificErrorAnalysis):

f := (x, t) -> x*t:
g := (x, t) -> x^2*t:

val := [seq(i,i=0.125..0.875,0.25)];

n     := numelems(val):
val_x := [seq(op(val[[k$n]]), k=1..n)]:
val_t := op~([val$n]):
val_f := f~(val_x, val_t):
val_g := g~(val_x, val_t):

res  := convert([val_t, val_f, val_g], Matrix)^+;

# Column titles as strings
cols := Vector[row](4, ["x", "t", "f(x, t)", "g(x, t)"]);

# Column titles as math symbols
# cols := Vector[row](4, ['x', 't', ''f(x, t)'', ''g(x, t)'']);

val := [.125, .375, .625, .875]

 

res := Vector(4, {(1) = ` 16 x 3 `*Matrix, (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

 

Vector[row](%id = 18446744078165679870)

(1)

TDI := style=TwoDimOutput:
FC  := fillcolor="LightGray":
W := Table(
           seq(Column(), j=1..n), widthmode=percentage, width=50, alignment=center,

           # Column titles as strings
           Row( seq( Cell(cols[j], FC), j=1..4) ),

           # Column titles as math symbols
           # Row( seq( Cell( Textfield(TDI, Equation(cols[j])), FC), j=1..4) ),

           seq(
             op([
               Row(
                    Cell( Textfield(TDI, Equation(val[k])), rowspan=4, padding=50),
                    seq(
                      Cell( Textfield(TDI, Equation(res[1+n*(k-1), j])) ),
                      j=1..3
                    )
               ),
               seq(
                 Row(
                      seq(
                        Cell( Textfield(TDI, Equation(res[i, j])) ),
                        j=1..3
                      )
                 ),
                 i=2+n*(k-1)..n+n*(k-1)
               )
             ]),
             k=1..n
           )
     ):
DocumentTools:-InsertContent(Worksheet(Group(Input( W )))):

 


 

Download Layout_1.mw

If you want to Export your data use the package ExcelTools (procedure Export).

It's still possible to save as a text file (package FileTools for instance) and next open this file from Excel by carefully choosing the forat of the data

You can also creat a spreadsheet within your MAPLE worksheet by using the Spread package.

Your code doesn' work with my 2015 version and OI have no time to adapt it.
So this is what you must do for this notional example
 

combs := combinat:-permute([a, b, c, d, e], 3)

arrangements := map(convert, {combs[]}, set)

 

First thing, a solution which uses all your code but the call to LeastSquaresPlot

restart:
UseHardwareFloats := false:

# your code

P := convert(points, listlist):
P := convert(P, Matrix):
Statistics:-LinearFit([1, x], P, x);
   1577.0067814989438 + 245.30627827585107 x


Second thing: your code is amazing slow to run;
You can improve Part A by observing that "dist" can be obtained by sampling a Binomial(n, 1/2) distribution

restart:
with(Statistics):

UseHardwareFloats := false:
dist := proc (n::nonnegint) 
  Sample(Binomial(n, 0.5), 1)[1]:
  round(%-(n-%))
end proc:
CodeTools:-Usage(dist(10^6))  # compare to your own "dist"


ave_dist is the mean of N (N=10000) realizations od dist(n).
A classical result in Statistics is this one :

  • if B is Binomial(n, p) and B' is Binomial(n', p) and B and B' are mutually independant, then 
    B+B' is Binomial(n+n', p)

Applied here this gives simply

ave_dist := proc (n::nonnegint, N::nonnegint) 
local B:
B := RandomVariable(Binomial(n*N, 0.5));
# what you want is (Sample(B, 1)[1] - Mean(B))/N
# or more simply
return ( Sample(B, 1)[1] - n*N/2 )/N
end proc:

And there is not need of procedure "dist" any longer (only ave_dist matters)
Now, let's build A:

n := 1000:
N := 10000:
M := 10:
A := Vector(M, i -> ave_dist(n*i, N))

Here is the complete code
 

restart;

with(Statistics):

UseHardwareFloats := true:

randomize():

ave_dist := proc (n::nonnegint, N::nonnegint)
local B:
B := RandomVariable(Binomial(n*N, 0.5));
# what you want is: ( Sample(B, 1)[1] - Mean(B) )/N
# or more simply
(Sample(B, 1)[1] - n*N/2)/N
end proc:

n := 1000:
N := 10000:
M := 10:
A := Vector(M, i -> ave_dist(n*i, N)):
 

V := Vector(M-1, i -> evalf(log(i+1))):

points := < <[0, 0]>^+, < V | V*~A[2..-1] > >:

fit := Statistics:-LinearFit([1, x], points, x);

-HFloat(0.00796016618590964)+HFloat(0.09968098510152414)*x

(1)

plots:-display(
  plot(fit, x=0..max(points[..,1]), color=blue, legend=typeset(evalf[4](fit))),
  plot(points, style=point, symbol=circle, symbolsize=15, legend="points"),
  size=[400, 400]
);

 

 


 

Download Faster.mw



The linear fit is of course very poor ... you are just trying to fit a random walk by a libear model...

 

Hi, 

This is not a cpomplete solution (see further) but an explanation: the problem comes from the equation eq2.
Considering it at the "theta equation" we expect (at least MAPLE does) fot it to be of the form 
theta''= F(theta', theta, + the other unknowns and their derivatives, eta) 
But it is fact of the form theta'' (a+b theta'') = F(...) which is not (at least in the classical theory of ODEs, but I do not know much from this one) an ODE.
Try just this to convince yourself

restart:
dsolve({ diff(x(t), t$2)*(1+diff(x(t), t$2))=diff(x(t), t)+1, x(0)=0, D(x)(0)=0}, numeric)
Error, (in DEtools/convertsys) unable to convert to an explicit first-order system

Again this is just an explanation.

But this simple example can be solved if you help MAPLE:
 

restart:
dsolve({ diff(x(t), t$2)*(1+diff(x(t), t$2))=diff(x(t), t)+1, x(0)=0, D(x)(0)=0}, numeric)
Error, (in DEtools/convertsys) unable to convert to an explicit first-order system

# isolate the term x''(t)
isolate(diff(x(t), t$2)*(1+diff(x(t), t$2))=diff(x(t), t)+1, diff(x(t), t$2));
                                                   (1/2)
         d  / d      \     1   1 /      / d      \\     
        --- |--- x(t)| = - - - - |5 + 4 |--- x(t)||     
         dt \ dt     /     2   2 \      \ dt     //     

# solve this new equation
dsolve({ %, x(0)=0, D(x)(0)=0}, numeric)

# of course it's possible that the solution becomes complex beyond some value 
# of t... but this is another problem

Explanation.mw

Applying this "isolation" strategy to your 4 equations makes the problem solvable.
Unfortunately an error occurs, probably because the discriminant the "isolated" eq2 contains is negative: 

 

restart; with(plots)

eq1 := 2*n*4^n*eta^((n+1)*(1/2))*(diff(f(eta), `$`(eta, 2)))^(n-1)*(diff(f(eta), `$`(eta, 3)))+4^n*(n+1)*eta^((n-1)*(1/2))*(diff(f(eta), `$`(eta, 2)))^n+4*f(eta)*(diff(f(eta), `$`(eta, 2)))-4*m*(diff(f(eta), eta))^2+m-2*M*(diff(f(eta), eta)) = 0; -1; eq1 := isolate(eq1, diff(f(eta), `$`(eta, 3)))

diff(diff(diff(f(eta), eta), eta), eta) = (1/2)*(-4^n*(n+1)*eta^((1/2)*n-1/2)*(diff(diff(f(eta), eta), eta))^n-4*f(eta)*(diff(diff(f(eta), eta), eta))+4*m*(diff(f(eta), eta))^2-m+2*M*(diff(f(eta), eta)))/(n*4^n*eta^((1/2)*n+1/2)*(diff(diff(f(eta), eta), eta))^(n-1))

(1)

eq2 := 2*eta*(diff(theta(eta), `$`(eta, 2)))+2*(diff(theta(eta), eta))+Pr*(f(eta)*(diff(theta(eta), eta))-s*(diff(f(eta), eta))*theta(eta))+Nb*(diff(theta(eta), eta))*(diff(phi(eta), eta))+Nt*(diff(theta(eta), `$`(eta, 2)))^2 = 0; -1; eq2 := isolate(eq2, diff(theta(eta), `$`(eta, 2)))

diff(diff(theta(eta), eta), eta) = (-eta+(Nt*theta(eta)*Pr*s*(diff(f(eta), eta))-Nb*Nt*(diff(phi(eta), eta))*(diff(theta(eta), eta))-Nt*f(eta)*Pr*(diff(theta(eta), eta))-2*Nt*(diff(theta(eta), eta))+eta^2)^(1/2))/Nt

(2)

eq3 := 2*eta*(diff(phi(eta), `$`(eta, 2)))+2*(diff(phi(eta), eta))+Sc*(f(eta)*(diff(phi(eta), eta))-s*(diff(f(eta), eta))*phi(eta))+Nb*(2*eta*(diff(theta(eta), `$`(eta, 2)))+2*(diff(theta(eta), eta)))/Nt = 0; -1; eq3 := isolate(eq3, diff(phi(eta), `$`(eta, 2)))
NULL

diff(diff(phi(eta), eta), eta) = (1/2)*(-2*(diff(phi(eta), eta))-Sc*(f(eta)*(diff(phi(eta), eta))-s*(diff(f(eta), eta))*phi(eta))-Nb*(2*eta*(diff(diff(theta(eta), eta), eta))+2*(diff(theta(eta), eta)))/Nt)/eta

(3)

eq4 := 2*eta*(diff(chi(eta), `$`(eta, 2)))+2*(diff(chi(eta), eta))+Lb*(f(eta)*(diff(chi(eta), eta))-s*(diff(f(eta), eta))*chi(eta))-Pe*(2*eta*chi(eta)*(diff(phi(eta), `$`(eta, 2)))+2*chi(eta)*(diff(phi(eta), eta))+2*eta*(diff(chi(eta), `$`(eta, 2)))*(diff(phi(eta), `$`(eta, 2)))) = 0; -1; eq4 := isolate(eq4, diff(chi(eta), `$`(eta, 2)))

diff(diff(chi(eta), eta), eta) = (-2*(diff(chi(eta), eta))-Lb*(f(eta)*(diff(chi(eta), eta))-s*(diff(f(eta), eta))*chi(eta))+2*chi(eta)*Pe*eta*(diff(diff(phi(eta), eta), eta))+2*chi(eta)*Pe*(diff(phi(eta), eta)))/(-2*Pe*eta*(diff(diff(phi(eta), eta), eta))+2*eta)

(4)

bcs := (D(f))(a) = 0, f(a) = 2*s*a*(D(phi))(a)/Sc, theta(a) = 1, phi(a) = 1, chi(a) = 1, (D(f))(10) = 1/2, theta(10) = 0, phi(10) = 0, chi(10) = 0:

params := {Lb = .1, M = .1, Nb = .6, Nt = .2, Pe = 5, Pr = 6.2, Sc = .1, a = 0.1e-1, m = 1/3, n = 1, s = .1};

{Lb = .1, M = .1, Nb = .6, Nt = .2, Pe = 5, Pr = 6.2, Sc = .1, a = 0.1e-1, m = 1/3, n = 1, s = .1}

(5)

print~(eval([eq1, eq2, eq3, eq4, bcs], params)):

diff(diff(diff(f(eta), eta), eta), eta) = (1/8)*(-8*(diff(diff(f(eta), eta), eta))-4*f(eta)*(diff(diff(f(eta), eta), eta))+(4/3)*(diff(f(eta), eta))^2-1/3+.2*(diff(f(eta), eta)))/eta

 

diff(diff(theta(eta), eta), eta) = -5.000000000*eta+5.000000000*(.124*(diff(f(eta), eta))*theta(eta)-.12*(diff(phi(eta), eta))*(diff(theta(eta), eta))-1.24*f(eta)*(diff(theta(eta), eta))-.4*(diff(theta(eta), eta))+eta^2)^(1/2)

 

diff(diff(phi(eta), eta), eta) = (1/2)*(-2*(diff(phi(eta), eta))-.1*f(eta)*(diff(phi(eta), eta))+0.1e-1*(diff(f(eta), eta))*phi(eta)-6.000000000*eta*(diff(diff(theta(eta), eta), eta))-6.000000000*(diff(theta(eta), eta)))/eta

 

diff(diff(chi(eta), eta), eta) = (-2*(diff(chi(eta), eta))-.1*f(eta)*(diff(chi(eta), eta))+0.1e-1*(diff(f(eta), eta))*chi(eta)+10*eta*chi(eta)*(diff(diff(phi(eta), eta), eta))+10*chi(eta)*(diff(phi(eta), eta)))/(-10*eta*(diff(diff(phi(eta), eta), eta))+2*eta)

 

(D(f))(0.1e-1) = 0

 

f(0.1e-1) = 0.2000000000e-1*(D(phi))(0.1e-1)

 

theta(0.1e-1) = 1

 

phi(0.1e-1) = 1

 

chi(0.1e-1) = 1

 

(D(f))(10) = 1/2

 

theta(10) = 0

 

phi(10) = 0

 

chi(10) = 0

(6)

infolevel[dsolve] := 4:

sol := dsolve(eval([eq1, eq2, eq3, eq4, bcs], params), numeric, output = listprocedure, maxmesh = 1024)

dsolve/numeric: entering dsolve/numeric

Error, (in dsolve/numeric/BVPSolve) unable to store '-HFloat(19.384294902530286)+HFloat(26.40565130500988)*I' when datatype=float[8]

 

``


 

Download Pblm2_isolated.mw




PS: these equations look like some boundary layer equations in fluid mechanics? If they are classical, algorithms do exist to solve them, have you done some research about this?

First 35 36 37 38 39 40 41 Last Page 37 of 52