Advanced Search

mmcdara

7891 Reputation

22 Badges

9 years, 55 days
Contact mmcdara

MaplePrimes Activity


These are answers submitted by mmcdara

Some typos in your pseudo code... here c...

mmcdara 7891
December 06 2020
1 2

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 manual...

mmcdara 7891
December 01 2020
0 0

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...

mmcdara 7891
December 01 2020
1 2

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:   ...

mmcdara 7891
November 30 2020
0 0

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 dou...

mmcdara 7891
November 29 2020
1 7

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 wit...

mmcdara 7891
November 29 2020
3 5


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 ...

mmcdara 7891
November 27 2020
0 0

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 v...

mmcdara 7891
November 27 2020
0 0

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 y...

mmcdara 7891
November 27 2020
3 0

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 solution but an e...

mmcdara 7891
November 26 2020
0 2

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?

Do you mean that you have two curves C1 ...

mmcdara 7891
November 26 2020
1 0

Do you mean that

  • you have two curves C1 and C2,
  • only known though a collection of N1 points for C1 and a collection N2 points for C2 
    (where N1 doesn't necessarily equals N2), 
  • and that you want to compute some distance between a function F1 whose graph would be C1 and a function F2 whose graph would be F2?

What I say here is something very close to what Acer said too.

If it's not that, there's no need to go further than the solution given by Carl.

Here is an illustration in the spirit of what acer said 

 

> 

restart:
> 

f := x -> cos(x):
g := x -> sin(x):
> 

NX := 10:
X  := Vector(NX, i -> 2.*Pi*(i-1)/(NX-1)):

NY := 8:
Y  := Vector(NY, i -> 2.*Pi*i/(NY-1)):
> 

F := < X | f~(X) >:
G := < Y | g~(Y) >:

F, G

Matrix(10, 2, {(1, 1) = 0., (1, 2) = 1., (2, 1) = .6981317009, (2, 2) = .7660444431, (3, 1) = 1.396263402, (3, 2) = .1736481773, (4, 1) = 2.094395103, (4, 2) = -.5000000005, (5, 1) = 2.792526804, (5, 2) = -.9396926211, (6, 1) = 3.490658504, (6, 2) = -.9396926208, (7, 1) = 4.188790205, (7, 2) = -.4999999998, (8, 1) = 4.886921906, (8, 2) = .1736481781, (9, 1) = 5.585053607, (9, 2) = .7660444435, (10, 1) = 6.283185308, (10, 2) = 1.}), Matrix(8, 2, {(1, 1) = .8975979011, (1, 2) = .7818314825, (2, 1) = 1.795195802, (2, 2) = .9749279122, (3, 1) = 2.692793703, (3, 2) = .4338837392, (4, 1) = 3.590391605, (4, 2) = -.4338837399, (5, 1) = 4.487989506, (5, 2) = -.9749279124, (6, 1) = 5.385587407, (6, 2) = -.7818314819, (7, 1) = 6.283185308, (7, 2) = 0.8204135232e-9, (8, 1) = 7.180783209, (8, 2) = .7818314830})

 (1)
> 

plots:-display(
  plot(convert(F, listlist), style=point, symbol=solidcircle, symbolsize=15, color=blue, legend="f"),
  plot(convert(G, listlist), style=point, symbol=solidcircle, symbolsize=15, color=red , legend="g"),
  plot(f(x), x=F[1, 1]..F[NX, 1], color=gray),
  plot(g(x), x=G[1, 1]..G[NY, 1], color=gray)
);

 
> 

Common_Range := max(X[1], Y[1])..min(X[-1], Y[-1])

.8975979011 .. 6.283185308

 (2)
> 

XP := [
        op(1, Common_Range),
        map(n -> if verify(X[n], Common_Range, 'interval') then X[n] end if, [$1..NX])[],
        op(2, Common_Range)
      ]:

YP := [
        op(1, Common_Range),
        map(n -> if verify(Y[n], Common_Range, 'interval') then Y[n] end if, [$1..NY])[],
        op(2, Common_Range)
      ]:
> 

#Acer's idea

SG := unapply( CurveFitting:-Spline(G, x, degree=1), x):
SF := unapply( CurveFitting:-Spline(F, x, degree=1), x):

#Carl's "maximum distance"
plot(abs(SF(x)-SG(x)), x=Common_Range, gridlines=true);
res     := maximize(abs(SF(x)-SG(x)), x=Common_Range, location=true):
DistMax := res[1];
LocMax  := eval(x, res[2][1][1]);


plots:-display(
  plot(convert(F, listlist), style=point, symbol=solidcircle, symbolsize=15, color=blue, legend="f"),
  plot(convert(G, listlist), style=point, symbol=solidcircle, symbolsize=15, color=red , legend="g"),
  plot(f(x), x=F[1, 1]..F[NX, 1], color=gray),
  plot(g(x), x=G[1, 1]..G[NY, 1], color=gray),
  plot([SG(x), SF(x)], x=Common_Range, gridlines=true, color=[red, blue]),
  plottools:-polygon(
    [ seq([x, SF(x)], x in XP),  seq([x, SG(x)], x in ListTools:-Reverse(YP)), [XP[1], SF(XP[1])] ],
    style='polygon', color=gold, transparency=0.7
  ),
  plot([[LocMax, SF(LocMax)], [LocMax, SG(LocMax)]], thickness=3, color=black)
);

 

1.378619853

  

5.385587407

  

 
> 

 


 

Download Maximum_Distance.mw

Change your last line by this one  ...

mmcdara 7891
November 25 2020
3 5

Change your last line by this one 
 

sol := dsolve(eval([ODEs[], bcs], params), numeric, output = listprocedure);

Error, (in dsolve/numeric/bvp/convertsys) too few boundary conditions: expected 8, got 7

And introduce a eighth boundary equation

Hi, Expressing "max" in t...

mmcdara 7891
November 25 2020
2 0

@Carl Love @Kitonum @vv

Hi, 

Expressing "max" in terms of Heaviside functions will give the result whether the graphs of f(x,p) and h(x,p) intersect or not:

> 

restart:
> 

f := (x, p) -> x + p;
g := (x, q) -> x^2 + q*x;

proc (x, p) options operator, arrow; x+p end proc

  

proc (x, q) options operator, arrow; x^2+q*x end proc

 (1)
> 

# h = max(f, g)

h := (x, p, q) -> f(x, p)*Heaviside(f(x, p)-g(x, q))+g(x, q)*Heaviside(g(x, q)-f(x, p))

proc (x, p, q) options operator, arrow; f(x, p)*Heaviside(f(x, p)-g(x, q))+g(x, q)*Heaviside(g(x, q)-f(x, p)) end proc

 (2)
> 

# example 1 (Katatonia's)

pq := [p=-1, q=1];

plot(eval([f(x, p), g(x, q), h(x, p, q)], pq), x=-2..4, color=[red, green, blue], legend=['f', 'g', 'h'])

[p = -1, q = 1]

  

 
> 

# a lengthy output not represented

hi := int(h(x, p, q), x=a..b):

# integration for Katatonia's example

hi := int(eval(h(x, p, q), pq), x=a..b);

(1/3)*b^3+(1/2)*b^2-(1/3)*a^3-(1/2)*a^2

 (3)
> 

# example continued

ab := [a=-2, b=3];

H := eval(h(x, op(rhs~(pq))), ab);

plots:-display(
  plot(H, x=-3..4, color=blue),
  plot(H, x=eval(a..b, ab), color=blue, filled=true),
  title=typeset('Int(h, a..b)'=eval(hi, ab))
);
  

[a = -2, b = 3]

  

(x-1)*Heaviside(-x^2-1)+(x^2+x)*Heaviside(x^2+1)

  

 
> 

# example 2

pq := [p=1, q=-1];

plot(eval([f(x, p), g(x, q), h(x, p, q)], pq), x=-2..4, color=[red, green, blue], legend=['f', 'g', 'h'])

[p = 1, q = -1]

  

 
> 

ab := [a=-1, b=3];

H := eval(h(x, op(rhs~(pq))), ab);

plots:-display(
  plot(H, x=-3..4, color=blue),
  plot(H, x=eval(a..b, ab), color=blue, filled=true),
  title=typeset('Int(h, a..b)'=eval(hi, ab))
);

[a = -1, b = 3]

  

(x+1)*Heaviside(-x^2+2*x+1)+(x^2-x)*Heaviside(x^2-2*x-1)

  

 
> 

 


 

Download With_Heaviside_instead_of_max.mw

Is it this that you want? (obtained wit...

mmcdara 7891
November 23 2020
1 0

Is it this that you want?

(obtained with MAPLE 2015, probably prettier with newer versions) 
 

> 

restart
> 

with(GraphTheory):
> 

with(SpecialGraphs):
> 

P := PetersenGraph():
DrawGraph(P);

 
> 

NC := ChromaticNumber(P, 'col');
V  := col;
MyColors := [ red, green, yellow];

3

  

[[1, 3, 8, 10], [2, 4, 6], [5, 7, 9]]

  

[red, green, yellow]

 (1)
> 

for n from 1 to NC do
  HighlightVertex(P, V[n], MyColors[n])
end do:
DrawGraph(P);

 
> 

 


 

Download colors.mw

Your problem is not precisely defined: W...

mmcdara 7891
November 22 2020
1 5

Your problem is not precisely defined:

  • What is the distribution of each element of the vector?
  • Are they identical distributions;
  • are they independant
  • what are they: Uniform, Gaussian...something else?


Depending on what you want, Kitonum may have given you the right answer.

Here is another reply
A classical problem in Statistics is the following:
"How to generate samples of K independant random variables X1, ..., XK, all with supports [0, 1], such that the sum of each P-uple is 1?"
This problem arises, for example, in the generation of mixtures of P chemical components.

The solution is given by the Dirichlet distribution  Dirichlet_distribution

Here is the code. I generate N=10 row vectors whose the sum of their elements is 1.
The parameter vector alpha is <1, 1, 1, 1>, which means that each component of a row vector has a (marginal) Uniform distribution with support [0, 1]; other choices of alpha correspond to other marginals.
Th
 

> 

restart:
> 

with(Statistics):
> 

#choose a vector of 4 strictly positive values

alpha := [1$4]:
> 

# sample 4 Beta random variables Gamma(1, alpha[k])

N := 10:
G := NULL:
for k from 1 to 4 do
  G||k := Sample(GammaDistribution(1, alpha[k]), N):
  G := G + G||k
end do:
> 

# compute N row vectors whose sum is 1

for k from 1 to 4 do
  X||k := G||k /~ G
end do:
> 

# verify

X := < seq(X||k, k=1..4) >:
< < [seq(cat("X", k), k=1..4), "sum"] >^+, < X, 4*~Mean(X) >^+ >

Matrix([["X1", "X2", "X3", "X4", "sum"], [0.184241082374523e-1, 0.711499184152927e-1, .162474048555286, .747951924791969, 1.], [.259867343263675, .300364851825076, .154427851679308, .285339953231941, 1.00000000000000], [.485558196598768, .234029799888272, .121682320144697, .158729683368264, 1.00000000000000], [.559789392921627, .286449768936391, 0.488868623759014e-1, .104873975766080, 1.], [0.857122439391480e-1, .414908231391176, .156425854494076, .342953670175600, 1.]])

 (1)
> 

 


 

Download DirichletDistribution.mw

 

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