Advanced Search

mmcdara

6568 Reputation

18 Badges

8 years, 103 days
Contact mmcdara

MaplePrimes Activity


These are answers submitted by mmcdara

Statistics:-LinearFit can do t...

mmcdara 6568
February 17 2019
1 0

Statistics:-LinearFit can do the job quite well. So there is no bug in Statistics:-LinearFit  (or at least not on this problem).

The issue comes from the fact that the problem is not relevant, per se, of the framework of Ordinary Least Squares (OLS) regression where we aim to assess the values of the parameters A and B of the model:

" Conditional expectation of Y given X=x = B*x+A "

The reason for this is that the slope is known.
In the attached file I show how a rewritting of the initial problem sets it back in the OLS framework. But this is neither immediate not natural.

However, the simplest method to assess the value of "a" comes from classical estimation results in OLS : the best unbiased estimator of "a" is just the mean of pts2 - 3*pts1

 


 

> 

restart:
> 

with(Statistics):
> 

pts1 := Vector([3, 5, 21]);
> 

pts2 := Vector([4, 6, 16])

pts1 := Vector(3, {(1) = 3, (2) = 5, (3) = 21})

  

pts2 := Vector(3, {(1) = 4, (2) = 6, (3) = 16})

 (1)

It couldn't be easier

> 

# From standard results in Ordinary Least Squares (OLS) regression,
# the Best Linear Unbiased Estimator (BLUE) of the intercept A
# in the model: Expectation(Y | X) = B*X+A is just
# BLUE(A) = Mean(Y) - BLUE(B)*Mean(X)
#
# Here BLUE(B) is supposed to be known and set to 3.
# Then

BLUE__A := Mean(pts2 - 3 *~ pts1);
add (pts2 - 3 *~ pts1) / 3

HFloat(-20.333333333333332)

  

-61/3

 (2)

Using Statistics:-LinearFit

> 

# The problem is not relevant of OLS but can be solved with LinearFit
# (at the price of a few pirouettes)
#
# Rewrite the problem this way:
#   1/ define Z as Y-3*X

Z := (X, Y) -> Y - 3*X;

#   2/ imagine you repeat 3 times the same experiment for a same value of the
#      control variable "a" ("a" is the only regressor of the problem)
#      Write Newpts2 the vector of the outcomes of these experiments
Newpts2 := Z~(pts1, pts2);
#   3/ define then the "design matrix" as an arbitrary constant vector,
#      for instance < 1, 1, 1>

Newpts1 := Vector([1, 1, 1]);

#   4/ apply LinearFit as usual, remembering that the only regressor is "a"

L := LinearFit([a], Newpts1, Newpts2, a);

#   5/ Substitute "a" by the value it takes in the vector Newpts1

A := subs(a=1, L);
convert(%, rational);

Z := proc (X, Y) options operator, arrow; Y-3*X end proc

  

Newpts2 := Vector(3, {(1) = -5, (2) = -9, (3) = -47})

  

Newpts1 := Vector(3, {(1) = 1, (2) = 1, (3) = 1})

  

-HFloat(20.33333333333334)*a

  

HFloat(-20.33333333333334)

  

-61/3

 (3)
> 

 

> 

 


 

Download Simplest.mw

Maybe this?...

mmcdara 6568
February 17 2019
1 3

Simulation of a "true" Galton board is very difficult.
You have, in all rigor, to simulate the drop of a sphere (a disc in 2D) shocking many fixed cylindrical pins (discs in 2D).
The trajectory of a disc (assuming the problem is 2D) is described by a succession of free flights submitted to the gravitational field (rather simple to code), each of them ending when the falling disc hurts a fixed pin (another disc).
Determining the impact point is a difficult problem.
More of this you have to describe the physics of the shock (elastic or not), account for possible slipping and rolling...
Last but not least a true Galton board has a limited width an the pins are confined in a kind of funnel whose walls modify the fall at at boundaries of the board.

What I propose you is a simpler simulation based on a random walk on a "triangular lattice" (Pascal's lattice).
Starting from the position (i, j) a "ball" can move either to position (i+1, j+1) or (i+1, j-1) with equal probabilities. The process is repeated from the initial position (0, 0) a given number of times.
Think to this to some kind of "discrete Galton board".

Look to the animation in the attached file.
If it suits you I will give you some explanation about it... if required


WATCH OUT: I faced difficulties to load the file with the animations, so I commented the corresponding line (red written on yellow background).
To run the animation, once the plot after the command Fallouts(100, 10) has been displayed, just clik on it and run the animation by using the keys in the upper left of your Maple worksheet.

It's likely that this code can be written in clever and more performing way.
 

> 

restart:

see for instance https://en.wikipedia.org/wiki/Random_walk

> 

with(Statistics):
with(plots):
with(plottools):
> 

randomize():
> 

R := rand(0 .. 1):

Fallouts := proc(NP, N)
   local  Fallout:
   local  GaltonBoard, depart, t, n, c, T, GBH;
   global cs, GB:
   
 
   # fallout of a single "ball"
 

   Fallout := proc(np)
      local GBH, depart, t, n, c, T;
      global cs, GB:


      GBH := copy(GB):
      depart := [0, 0];
      t      := depart:
      for n from 1 to N do
         depart := depart +~ [1, (2*R()-1)];
         t      := t, depart
      end do:   
      t  := [t, [depart[1]+1, depart[2]]];
      c  := unapply(CurveFitting:-Spline(op~(1, t), op~(2, t), x, degree=1), x);
      T  := transform((x,y) -> [y+N+1, x]):
  
      if np > 1 then
         GBH := display(GBH, T(Histogram([cs], discrete, thickness=5, frequencyscale=absolute)) )
      end if;
      
      cs := cs, depart[2];
   
   plots[animate](plot, [c(x), x=0..A, color=red, thickness=2], A=0..N+1, frames=N+2, background=GBH);
   
   end proc:

   GaltonBoard := [seq(seq([i, j], j in [seq(-i..i, 2)]), i=0..N)]:
   GB := plot(GaltonBoard, style=point, symbol=solidcircle, symbolsize=10, scaling=constrained, gridlines=true, color=blue):

   cs := NULL:

   plots[animate]( Fallout, [np], np=1..NP, frames=NP);

end proc:
> 

global cs:

# drop 100 "balls" on a "discrete Galton Board" made oh 10 layers of pins

# Fallouts(100, 10);   # TO BE UNCOMMENTED
> 

cs

cs

 (1)
> 

# Histogram([cs])
> 

 


 

Download Galton-RandomWalk.mw

Maybe you could try this   ...

mmcdara 6568
February 07 2019
2 5

Maybe you could try this


 

> 

restart
> 

with(Statistics):
> 

XG := h -> RandomVariable(Normal(0, h))

proc (h) options operator, arrow; Statistics:-RandomVariable(Normal(0, h)) end proc

 (1)
> 

MaxIt := 1000

1000

 (2)
> 

randomize():
> 

# Be careful: h must be strictly positive

h := evalf([seq(1-k/MaxIt, k=0..MaxIt-1)]):
h[1],h[-1];

1., 0.1000000000e-2

 (3)
> 

# do you really need to print the results ?
#
# for i from 1 to MaxIt do
#    printf("%a\n%a\n", Sample(XG(h[i]), 5), i);
# end do:
> 

# Maybe this could suit you?


# Stage 1 ; generate a random matrix by sampling Normal(0, 1)

M := Sample(XG(1), [MaxIt, 5], method = [envelope, range = -1 .. 1]):


# Watch out: stage 2 is computationnaly expensive.
#            It is presented here to help understanding the method
#
# Stage 2: rescale M according to the standard deviations 1, 0.999, ...0.001

S := LinearAlgebra:-DiagonalMatrix(h):

# Stage 3: do the product S.M

Res := S.M

Res := Vector(4, {(1) = ` 1000 x 5 `*Matrix, (2) = `Data Type: `*float[8], (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

 (4)
> 

# The same thing in a more efficient way

M   := CodeTools:-Usage( Sample(XG(1), [MaxIt, 5], method = [envelope, range = -1 .. 1]) ):
Res := CodeTools:-Usage( <  seq(M[n, ..] *~ h[n], n=1..MaxIt) > );

memory used=0.73MiB, alloc change=0 bytes, cpu time=35.00ms, real time=11.00ms, gc time=0ns
memory used=0.60MiB, alloc change=0 bytes, cpu time=9.00ms, real time=3.00ms, gc time=0ns

  

Res := Vector(4, {(1) = ` 1000 x 5 `*Matrix, (2) = `Data Type: `*float[8], (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

 (5)
> 

for i from 1 to MaxIt do
  # printf("%3d : %a\n", i, [entries(Res[i,..], nolist)]);
end do:
> 

 


 

Download Gaussian_Reply.mw

tomleslie has provide a clever and fast ...

mmcdara 6568
January 31 2019
0 4

Several people have already answered your question in different ways

To complete their propositions, here is a time and memory comparison of different ways to answer yout problem

> 

restart:
> 

with(Statistics):   # the most natural (?) and one of the fastest if you want 0 and 1 outputs
B := RandomVariable(Bernoulli(1/2)):
> 

numThrows:= 10^5:
> 

CodeTools:-Usage( Sample(B, numThrows) ):

memory used=0.87MiB, alloc change=0 bytes, cpu time=7.00ms, real time=7.00ms, gc time=0ns

  
> 

restart;
> 

with(StringTools) : # a simple way to get "H" and "T" ouputs
> 

numThrows:= 10^5:
CodeTools:-Usage( Random(numThrows, "HT" ) ):                     # a string of the form "HHT..."
CodeTools:-Usage( Explode(Random(numThrows, "HT" )) ):            # a list of the form [ "H", "H", "T", ...
CodeTools:-Usage( Vector(Explode(Random(numThrows, "HT" ))) ):    # a vector (here fot comparison with tomleslie's solution below)

memory used=154.93KiB, alloc change=0 bytes, cpu time=193.00ms, real time=193.00ms, gc time=0ns
memory used=0.86MiB, alloc change=0 bytes, cpu time=192.00ms, real time=193.00ms, gc time=0ns

memory used=22.27MiB, alloc change=0 bytes, cpu time=383.00ms, real time=384.00ms, gc time=41.38ms

  
> 

restart;
> 

# tomleslie  ... probably the fastest (?) to deliver "H" and "T" ouputs
numThrows:= 10^5:
> 

r:= rand(1..2):
> 

CodeTools:-Usage(  Vector[row](numThrows, i->r()) ):
CodeTools:-Usage(  Vector[row](numThrows, i->["H","T"][r()]) ):

memory used=0.81MiB, alloc change=0 bytes, cpu time=33.00ms, real time=33.00ms, gc time=0ns
memory used=0.77MiB, alloc change=0 bytes, cpu time=38.00ms, real time=38.00ms, gc time=0ns

  
> 

restart;
> 

with(LinearAlgebra): # advantage: there exist a very fast way to tansform 0 & 1 into "H" & "T"
> 

numThrows:= 10^5:
> 

V := CodeTools:-Usage(  RandomVector[row](numThrows,generator=rand(0..1)) ):  # outputs are 0 or 1
W := CodeTools:-Usage(  subs({0="H", 1="T"}, V) ):                            # outputs are "H" or "T"

memory used=4.70MiB, alloc change=2.00MiB, cpu time=50.00ms, real time=51.00ms, gc time=5.63ms
memory used=0.76MiB, alloc change=0 bytes, cpu time=3.00ms, real time=3.00ms, gc time=0ns

  
> 

restart;
> 

with(SignalProcessing):  # comparable to Statistics
> 

numThrows:= 10^5:
> 

V := CodeTools:-Usage( floor~(GenerateUniform(numThrows, 0, 2)) ):  # the fastest (?), use Hfloats (default)

UseHardwareFloats := false:
CodeTools:-Usage( floor~(GenerateUniform(numThrows, 0, 2)) ):  # don't use Hfloats

memory used=1.79MiB, alloc change=0 bytes, cpu time=7.00ms, real time=8.00ms, gc time=0ns

memory used=203.74MiB, alloc change=47.02MiB, cpu time=1.58s, real time=1.59s, gc time=88.87ms

  
> 

restart:
> 

with(combinat):  # just for fun
> 

numThrows:= 10^5:
> 

CodeTools:-Usage( randperm(["T"$numThrows, "H"$numThrows])[1..numThrows] ):

memory used=20.67MiB, alloc change=4.58MiB, cpu time=161.00ms, real time=162.00ms, gc time=61.65ms

  
> 

 


 

Download Still_Another_Way.mw

 

After some digressions unrelated to your...

mmcdara 6568
January 28 2019
0 0

After some digressions unrelated to your initial question, here is my definitive contribution

final.mw

> 

restart:
> 

with(Statistics):
> 

Seed := randomize():
> 

alias(RV=RandomVariable):
> 

# Number of dice

N := 3;

3

 (1)
> 

# These two lines can be replaced by any algorithm amid those presented in this
# same thread (just choose the one you want)

S := add(RV(DiscreteUniform(1, 6)), k=1..N):
P := Probability~(S =~ [$N..6*N])

[1/216, 1/72, 1/36, 5/108, 5/72, 7/72, 25/216, 1/8, 1/8, 25/216, 7/72, 5/72, 5/108, 1/36, 1/72, 1/216]

 (2)
> 

# Define a discrete RV with outcomes of probabilities given by P
#
# Watch out: by default (I do not know if and how this could be changed), this RV has outcomes 1..5*N+1
# This raises an issue that will be fixed when the resuklts are ploted (the outcomes are N..6*N)

Dice := RV(ProbabilityTable(P))

_R2

 (3)
> 

# Sample "Dice" as any other RV

M := 100:
X := Sample(Dice, M):
> 

h := Histogram(X, minbins=5*N+1):

# The next operation is aimed to superimpose the exact mss function (aka P) and its
# empirical estimation h.
#
# I transform h into a list that can be plot with ColumnGraph.
# Watch out:
#      Maple 2018: use op~(2, op~(1, [op(1..-2, op(1, h))] ));
#      Maple 2015: use map(u -> print(u[-1][2]), [op(1..-2, op(1, h))]);

interface(version);

H := map(u -> u[-1][2], [op(1..-2, op(1, h))]):

ColumnGraph(
    [H, P],
    legend=["Sampled", "Exact"],
    color=[red, yellow],
    background=gray,
    title=cat("Sample size = ", M),
    tickmarks=[[seq(i=i+N, i=0..5*N)], default],
    size=[600, 400]
)

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

  

 
> 

 


 

Download final.mw

 

 

@tomleslie  A few remarks: Use&...

mmcdara 6568
January 26 2019
1 1

@tomleslie 

A few remarks:

  1. Use Quantile(X, k+i+j, numeric) for a faster generation:
    CodeTools:-Usage(Matrix(100, 10, (i,j) -> Quantile(X, (i-1)/100+(j-1)/1000, numeric))):
    memory used=14.25MiB, alloc change=0 bytes, cpu time=158.00ms, real time=159.00ms, gc time=0ns
    CodeTools:-Usage(Matrix(100, 10, (i,j) -> Quantile(X, (i-1)/100+(j-1)/1000))):
    memory used=430.05MiB, alloc change=0 bytes, cpu time=6.28s, real time=5.90s, gc time=555.98ms

    Note the CPU time (in this particular case) can be halved: Quantile(X, p) = -Quantile(X, 1-p)
     
  2. Avoid computing Quantile(X, 0) which is obviously -infinity and not -15.681...
     

Here is a version which uses DocumentTools for a smarter rendering

QuantileTable.mw


 

As you wrote it, your maplet can...

mmcdara 6568
January 05 2019
0 1

As you wrote it, your maplet can't work: Evaluate needs a function as argument, not a sequence of instructions.
Here is the correct way to code it in a simple case (no "for a in .... end do" loop ... but you will be able to generalize with from Carl's answer).

PS: the code I wrote conforms to your's, but it could be improved (remark first that there is no need to click on the "Display function" function button to visualize the MathML code in the MMLV field). It's also better to define "Input function  f(x)=" within a Label component, ...

> 

restart:
> 

with(Maplets[Elements]):
with(Student[Calculus1]):
> 

f := proc()
local g, TL2:

g := proc(func, slp)
  local a;
  a := Tangent(func,x=slp):
  Maplets:-Tools:-Set('MMLV1'=a)
end proc:


TL2:=Maplet([
  [
    "Input function  f(x)=", TextField['func'](20)
  ],
  [
    Button("Display function", Action(Evaluate('MMLV'='MathML[Export](func)')))
  ],
  MathMLViewer['MMLV'](),
  [
    "Input slope", TextField['slp'](10)
  ],  
  [
    Button("Tangent line", Action( Evaluate('function'=g, Argument(func), Argument(slp))))
  ],   
  MathMLViewer['MMLV1']()   
]):

Maplets[Display](TL2):

end proc:

f()
> 

 


 

Download Maplet.mw

 

your mistake...

mmcdara 6568
November 01 2018
1 0

You had forgotten that the projected dx*dy surface wasn't the surface itself:


 

> 

restart;
> 

with(VectorCalculus):
> 

SetCoordinates(cartesian[x, y, z]):

> 

R := sqrt(2):
s := x^2+y^2+z^2-R^2;
u := VectorField([x, 1, z]);

> 

N := Gradient(s);

x^2+y^2+z^2-2

  

u := Vector(3, {(1) = x, (2) = 1, (3) = z}, attributes = [vectorfield, coords = cartesian[x, y, z]])

  

N := Vector(3, {(1) = 2*x, (2) = 2*y, (3) = 2*z}, attributes = [vectorfield, coords = cartesian[x, y, z]])

 (1)
> 

n := N/sqrt(add(N[k]^2, k = 1 .. 3));
u . n;

un := algsubs(x^2+y^2+z^2=R^2, %):
un := simplify(un)

n := x*`#mover(mi("e"),mo("&lowbar;"))`[x]/sqrt(x^2+y^2+z^2)+y*`#mover(mi("e"),mo("&lowbar;"))`[y]/sqrt(x^2+y^2+z^2)+z*`#mover(mi("e"),mo("&lowbar;"))`[z]/sqrt(x^2+y^2+z^2)

  

x^2/(x^2+y^2+z^2)^(1/2)+y/(x^2+y^2+z^2)^(1/2)+z^2/(x^2+y^2+z^2)^(1/2)

  

-(1/2)*2^(1/2)*y^2+(1/2)*2^(1/2)*y+2^(1/2)

 (2)
> 

intr1 := solve(subs(z = 0, s), y);
> 

intr2 := solve(subs(z = 0, y = 0, s), x);

(-x^2+2)^(1/2), -(-x^2+2)^(1/2)

  

2^(1/2), -2^(1/2)

 (3)
> 

# forgot that the projected dx*dy surface is not the surface irself

l     := sqrt(x^2+y^2);
h     := sqrt(2-x^2-y^2);
phi   := arctan(h/l);
scale := simplify(1/sin(phi)) assuming x^2+y^2 > 0;

# surface of the sphere

R1 := 2*int(int(scale, y = intr1[2] .. intr1[1]), x = intr2[2] .. intr2[1], numeric);

evalf(4*Pi*R^2);

(x^2+y^2)^(1/2)

  

(-x^2-y^2+2)^(1/2)

  

arctan((-x^2-y^2+2)^(1/2)/(x^2+y^2)^(1/2))

  

2^(1/2)/(-x^2-y^2+2)^(1/2)

  

25.13274122

  

25.13274123

 (4)
> 

R1 := 2*Student:-MultivariateCalculus:-MultiInt(scale*un, y = intr1[2] .. intr1[1], x = intr2[2] .. intr2[1])

(16/3)*2^(1/2)*Pi

 (5)
> 

evalf(%);

23.69537567

 (6)
> 

du := Divergence(u)

2

 (7)
> 

intr1 := solve(s, z);
> 

intr2 := solve(subs(z = 0, s), y);
> 

intr3 := solve(subs(z = 0, y = 0, s));
> 

R2 := int(int(int(du, z = intr1[2] .. intr1[1]), y = intr2[2] .. intr2[1]), x = intr3[2] .. intr3[1]);
> 

evalf(%)

(-x^2-y^2+2)^(1/2), -(-x^2-y^2+2)^(1/2)

  

(-x^2+2)^(1/2), -(-x^2+2)^(1/2)

  

2^(1/2), -2^(1/2)

  

(16/3)*2^(1/2)*Pi

  

23.69537567

 (8)
> 

 


 

Download GaussXYZ.mw

In your particular case it's far mor...

mmcdara 6568
November 01 2018
1 0

In your particular case it's far more easier to verify the Gass theorem in spherical coordinates

I'll give a closer look to your code later
 

> 

restart;
> 

with(IntegrationTools):
with(VectorCalculus):
SetCoordinates(cartesian[x, y, z]):

> 

MySphere := x^2+y^2+z^2=2;
u := VectorField([x, 1, z]);

> 

N := Gradient(lhs(MySphere));

x^2+y^2+z^2 = 2

  

u := Vector(3, {(1) = x, (2) = 1, (3) = z}, attributes = [vectorfield, coords = cartesian[x, y, z]])

  

N := Vector(3, {(1) = 2*x, (2) = 2*y, (3) = 2*z}, attributes = [vectorfield, coords = cartesian[x, y, z]])

 (1)
> 

n := N/sqrt(add(N[k]^2, k = 1 .. 3));
n := subs(MySphere, n);

n := x*`#mover(mi("e"),mo("&lowbar;"))`[x]/sqrt(x^2+y^2+z^2)+y*`#mover(mi("e"),mo("&lowbar;"))`[y]/sqrt(x^2+y^2+z^2)+z*`#mover(mi("e"),mo("&lowbar;"))`[z]/sqrt(x^2+y^2+z^2)

  

n := Vector(3, {(1) = (1/2)*2^(1/2)*x, (2) = (1/2)*2^(1/2)*y, (3) = (1/2)*2^(1/2)*z}, attributes = [vectorfield, coords = cartesian[x, y, z]])

 (2)
> 

un := DotProduct(u, n);
du := Divergence(u);

(1/2)*x^2*2^(1/2)+(1/2)*2^(1/2)*y+(1/2)*z^2*2^(1/2)

  

2

 (3)
> 

SphereRadius := sqrt(rhs(MySphere));

2^(1/2)

 (4)
> 

J := (r, phi) -> r^2*sin(phi);

int(int(int(du*J(r, phi), theta=0..2*Pi), phi=0..Pi), r=0..SphereRadius);
evalf(%);

proc (r, phi) options operator, arrow; VectorCalculus:-`*`(r^2, sin(phi)) end proc

  

(16/3)*Pi*2^(1/2)

  

23.69537567

 (5)
> 

int(int(subs(r=SphereRadius, unsph)*J(SphereRadius, phi), theta=0..2*Pi), phi=0..Pi);
evalf(%);

(16/3)*Pi*2^(1/2)

 (6)
> 

 


 

Download Gauss.mwGauss.mw

Not the solution, just a way to derive t...

mmcdara 6568
October 27 2018
2 4

Not the solution, just a way to derive the equations by applying differential operators.
I think you have an ambiguity about the pressure (you declare p bar as a constant but it has a priori different values at z=-1 and z=+1 (unless pa is null)

For the solution: I'm not familiar with pdsolve and PDETools, but I'm afraid that you will have to do the things by hand. 
If I'm not mistaken this is a classical poiseuille flow with uniform velocity at the inlet (parabolic profile expected at the outlet) and I remember solving it by hand long ago at school. Maybe you will have to reproduce explicitely the same modus operandi in Maple?

Wait for a more enlightened reply about pdsolve/PDEtools


 

> 

restart:
> 

with(Student[VectorCalculus]):
> 

SetCoordinates(cylindrical[r,theta, z]):
> 

alias(u=u(r, 'theta', z));
alias(v=v(r, 'theta', z));
alias(p=p(r, 'theta', z));

u

  

u, v

  

u, v, p

 (1)
> 

V := VectorField(< v, 0, u >);
Divergence(V);

V := Vector(3, {(1) = v(r, theta, z), (2) = 0, (3) = u(r, theta, z)}, attributes = [vectorfield, coords = cylindrical[r, theta, z]])

  

(v+r*(diff(v, r))+r*(diff(u, z)))/r

 (2)
> 

Continuity_equation := eval(expand(%));

v/r+diff(v, r)+diff(u, z)

 (3)
> 

Laplacian(V):
Diffusion_term := beta^2 *~ map(expand, eval(%));

Diffusion_term := Vector(3, {(1) = beta^2*(-v(r, theta, z)/r^2+(diff(v(r, theta, z), r))/r+diff(diff(v(r, theta, z), r), r)+(diff(diff(v(r, theta, z), theta), theta))/r^2+diff(diff(v(r, theta, z), z), z)), (2) = 2*beta^2*(diff(v(r, theta, z), theta))/r^2, (3) = beta^2*(diff(diff(u(r, theta, z), z), z)+(diff(u(r, theta, z), r))/r+diff(diff(u(r, theta, z), r), r)+(diff(diff(u(r, theta, z), theta), theta))/r^2)})

 (4)
> 

Pressure_term := 1 *~ Gradient(p)

Pressure_term := Vector(3, {(1) = diff(p(r, theta, z), r), (2) = (diff(p(r, theta, z), theta))/r, (3) = diff(p(r, theta, z), z)})

 (5)
> 

Advection_term := 1 *~ Vector[column](3, i -> DotProduct(V, Gradient~(V[i])))

Advection_term := Vector(3, {(1) = v(r, theta, z)*(diff(v(r, theta, z), r))+u(r, theta, z)*(diff(v(r, theta, z), z)), (2) = 0, (3) = v(r, theta, z)*(diff(u(r, theta, z), r))+u(r, theta, z)*(diff(u(r, theta, z), z))})

 (6)
> 

C := :-Vector[column](3, i -> Advection_term[i] = - Pressure_term[i] + Diffusion_term[i])[[1, 3]]

C := Vector(2, {(1) = v(r, theta, z)*(diff(v(r, theta, z), r))+u(r, theta, z)*(diff(v(r, theta, z), z)) = -(diff(p(r, theta, z), r))+beta^2*(-v(r, theta, z)/r^2+(diff(v(r, theta, z), r))/r+diff(diff(v(r, theta, z), r), r)+(diff(diff(v(r, theta, z), theta), theta))/r^2+diff(diff(v(r, theta, z), z), z)), (2) = v(r, theta, z)*(diff(u(r, theta, z), r))+u(r, theta, z)*(diff(u(r, theta, z), z)) = -(diff(p(r, theta, z), z))+beta^2*(diff(diff(u(r, theta, z), z), z)+(diff(u(r, theta, z), r))/r+diff(diff(u(r, theta, z), r), r)+(diff(diff(u(r, theta, z), theta), theta))/r^2)})

 (7)
> 

NS := [Continuity_equation, entries(C, nolist)]

[v/r+diff(v, r)+diff(u, z), v*(diff(v, r))+u*(diff(v, z)) = -(diff(p, r))+beta^2*(-v/r^2+(diff(v, r))/r+diff(diff(v, r), r)+(diff(diff(v, theta), theta))/r^2+diff(diff(v, z), z)), v*(diff(u, r))+u*(diff(u, z)) = -(diff(p, z))+beta^2*(diff(diff(u, z), z)+(diff(u, r))/r+diff(diff(u, r), r)+(diff(diff(u, theta), theta))/r^2)]

 (8)
> 

# You can say "p bar is a constant and { p(z=-1)=0 & p(z=-1)=0 & p(z=-+1)=p__a } unless p__a=0 too  

BC := [
:-D[1](u)(0, 'theta', z) = 0,
v(0, 'theta', z) = 0,
u(1, 'theta', z) = 0,
v(1, 'theta', z) = p+alpha
]:

seq(print(BC[i]), i=1..numelems(BC))

(D[1](u))(0, theta, z) = 0

  

v(0, theta, z) = 0

  

u(1, theta, z) = 0

  

v(1, theta, z) = p+alpha

 (9)
> 

 


 

Download NavierStokes.mw

Replace your line   Diam; by  ...

mmcdara 6568
October 25 2018
0 1

Replace your line   Diam;
by                           Diam := simplify(Diam);

Nevertheless I would be tempted to say that Maple should simplify automatically.

Your line evalf(Diam) is corrupted.
Just rewrite it

TelescopeUnits_reply.mw


PS:Personnaly I never use this document mode input that I find very capricious and subjected to hide characters mistakenly typed

I'm not sure I understand you perfec...

mmcdara 6568
October 19 2018
0 1

I'm not sure I understand you perfectly.
What I think is: you have a discrete RV "X" with mass function fX(x)=1/(x+1) for x = 0, ..., k, where k is some given number (>0).

If It is that, a way to proceed is to use a ProbabilityTable Distribution

 

> 

restart:
> 

with(Statistics):
> 

# Discrete distribution from a probability table
# The sum of all masses must be equal to 1.

P := k -> [seq(1/(x+1), x=0..k)] /~ add(1/(x+1), x=0..k) ;

# example:

K := 5:
N := 10^3:

P(K);

X := RandomVariable(ProbabilityTable(P(K))):

S := Sample(X, N):
S[1..10];
Histogram(S);

 

P := proc (k) options operator, arrow; `~`[:-`/`]([seq(1/(x+1), x = 0 .. k)], ` $`, add(1/(x+1), x = 0 .. k)) end proc

  

"[[Typesetting:-mfrac(Typesetting:-mn("20",mathvariant = "normal"),Typesetting:-mn("49",mathvariant = "normal"),linethickness = "1",denomalign = "center",numalign = "center",bevelled = "false"), Typesetting:-mfrac(Typesetting:-mn("10",mathvariant = "normal"),Typesetting:-mn("49",mathvariant = "normal"),linethickness = "1",denomalign = "center",numalign = "center",bevelled = "false"), Typesetting:-mfrac(Typesetting:-mn("20",mathvariant = "normal"),Typesetting:-mn("147",mathvariant = "normal"),linethickness = "1",denomalign = "center",numalign = "center",bevelled = "false"), Typesetting:-mfrac(Typesetting:-mn("5",mathvariant = "normal"),Typesetting:-mn("49",mathvariant = "normal"),linethickness = "1",denomalign = "center",numalign = "center",bevelled = "false"), Typesetting:-mfrac(Typesetting:-mn("4",mathvariant = "normal"),Typesetting:-mn("49",mathvariant = "normal"),linethickness = "1",denomalign = "center",numalign = "center",bevelled = "false"), Typesetting:-mfrac(Typesetting:-mn("10",mathvariant = "normal"),Typesetting:-mn("147",mathvariant = "normal"),linethickness = "1",denomalign = "center",numalign = "center",bevelled = "false")]]"

  

Vector[row]([3., 2., 6., 2., 1., 1., 2., 1., 2., 1.])

  

 
> 

 


 

Download DiscreteRV.mw

Hi,  If your only concern is ...

mmcdara 6568
October 19 2018
2 3

Hi, 

If your only concern is to compare the model and the experiment, and if this experiment is a "not too large" collection of triples (t, x, y) , then the option output=Array[ ...], where ... is the list of the experimental t values, is probably the simpler solution.

Regarding now your "comparison metric".
The answer is very simple: there is unfortunately no good way to proceed!
Computing "the Root Mean Square of all (xexperiment(t)-xmodel(t), yexperiment(t)-ymodel(t))" seems very reasonable.
But it can be a very poor solution for some problems, specifically if x or y are of different natures, vary on very different scales and so on.
Luckily it doesn't seem to be the case here.

Another point you must be aware of: if the (t, x, y) experimental values come from measurements, there are probably entailed by "errors". Suppose that x and y are perfectly measured, but that t is not.
For instance the "true" t is T=1, and the measured value is t=1.1: then comparing the experiment at time t=1.1 to the model at time t=1.1 doesn't make sense.
A (the?) correct approach in this case is to use a bayesian framework where one infer the true value of T given its measure t.

Maybe you could be interested in "fitting" your model on the experiments? For instance to assess the value of A ...?
You also mention that "The problem is however that even if the (x,y) curves are rather close, if may not be a good model".
Are you concerned with aapossible "model error"?
Here again, for these two questions, the correct answer relies to the "bayesian calibration of computer model" issue.

Let me know if you are interested by taht.

To conclude: try your  "Root Mean Square of all (xexperiment(t)-xmodel(t), yexperiment(t)-ymodel(t)", examine carefully the residues (the model-to-experiment distance at each time t) and if all looks good, then it's done 
 

MAT := < header1, header2, <Column...

mmcdara 6568
October 18 2018
2 0

MAT := < header1, header2, <Column(t1[2,1],1..2) | Column(t2[2,1],2) > >:
interface(rtablesize=max(Dimensions(MAT))):
MAT;

 

 

Assuming the DEFINITION below suits you,...

mmcdara 6568
September 30 2018
0 3

Assuming the DEFINITION below suits you, you will find some elements in the attached file.

DEFINITION:

  • Let E the ellipse defined by the equation (x/a)^2+(y/b)^2=1
  • For each point P on E one define s(P) as the measure of the angle (counter clockwise counted) between Ox and OM
  • Let A and B two points on E such that
    • s(A) > =  0
    • s(B) > s(A)
    • s(B) <= 2*Pi
  • Let M any point such that s(A) <= s(M) <= s(B)
  • Let x(A) the abscissa of A and x(B) the one of B

Then the mean value of the distance d(O, M) from O=(0, 0) to M as M moves from A to B is defined by the integral of d(O, M) between x(A) and x(B), divided by x(B)-x(A).

In file Ellipse2.mw you will find a procedure 'md' which computes the mean value of d(O, M). It takes 4 arguments:

  1. a 
  2. b
  3. alpha = s(A)
  4. beta = s(B)

'md' handles (if I'm not mistaken) the 3 cases:

  1. alpha < Pi and beta < Pi
  2. alpha < Pi and beta > Pi
  3. alpha > Pi and beta > Pi 


Before that you will find a few lines that treat more formally the first situation alpha < Pi and beta < Pi

As claimed in this file, I was not capable to obtain a general formal solution for any values of a and b (that's why I arbitrarily wrote a:=2 and b:=1 at the top at the file).
The search of a general formal solution requires using assumptions about a, b, p, q (see the commands in the file) and the 'int' procedure of Maple is not very comfortable with 'assume' or 'assuming' statements (someone here will probably correct me on this point)

Ellipse2.mw

First 49 50 51 52 53 54 55 Page 51 of 56