Advanced Search

mmcdara

7891 Reputation

22 Badges

9 years, 55 days
Contact mmcdara

MaplePrimes Activity


These are answers submitted by mmcdara

@Michele95 As Carl Love already sai...

mmcdara 7891
November 22 2020
0 0

@Michele95 

As Carl Love already said it's difficult to find the location of the minima while N, h__TOT and t are not fixed (see the first part of the attached file).
But this is much simpler once thes parameters have been instanciated (look after the plot).
Here the solution is obtained by solving for h diff(alpha, h)=0 and checking the sign of the second derivative.

As your function may have discontinuities (for the triple (t, N, h__TOT) I chosed they are located at 0, 1, 5), minimize cannot find the solution unless you help it.
I wasn't capable to obtain a solution with Optimization:-Minimize neither (but maybe others could).


 

> 

restart:
> 

W1:= h-> piecewise(h < 1, 1, 2);
> 

W2:= h-> piecewise(h < 1, 3, 4);

proc (h) options operator, arrow; piecewise(h < 1, 1, 2) end proc

  

proc (h) options operator, arrow; piecewise(h < 1, 3, 4) end proc

 (1)
> 

alpha := (1/2*W1(h)*t+(W2(h)+N)*(t+1/2*t*h/(h__TOT-h)))/(h/2)/(W1(h)+W2(h))

alpha := (2*((1/2)*piecewise(h < 1, 1, 2)*t+(piecewise(h < 1, 3, 4)+N)*(t+(1/2)*t*h/(h__TOT-h))))/(h*(piecewise(h < 1, 1, 2)+piecewise(h < 1, 3, 4)))

 (2)
> 

# locations of extrema
#
# watchout: the values of t, N and h__TOT must guarantee that
# the two values in left_sol are < 1 and the two in right_sol are > 1


da := diff(alpha, h):
left_sol  := [ solve(da=0, h) ] assuming h < 1;
right_sol := [ solve(da=0, h) ] assuming h > 1;

[(2*N+7+(2*N^2+13*N+21)^(1/2))*h__TOT/(4+N), (2*N+7-(2*N^2+13*N+21)^(1/2))*h__TOT/(4+N)]

  

[(2*N+10+(2*N^2+18*N+40)^(1/2))*h__TOT/(N+6), (2*N+10-(2*N^2+18*N+40)^(1/2))*h__TOT/(N+6)]

 (3)
> 

# To decide if an extremal point gives a maximum or a minimum one must
# look the value of the second derivative of alpha at this extremal point.
# But one must also give numerical values to the remaining parameters,
# for instance

data := [N=10, h__TOT=5, t=3]:

d2a := diff(eval(da, data), h):
 
> 

# Firstly select the extremal points that are < than  1 or > than 1
LS := select(`<`, evalf(eval(left_sol , data)), 1);
RS := select(`>`, evalf(eval(right_sol , data)), 1);
 

[]

  

[15.77934423, 2.970655772]

 (4)
> 

# Next eval d2a at LS and RS

map(u -> evalf(eval(d2a, h=u)), RS);

[-0.3541865518e-2, 2.819541866]

 (5)
> 

# only RS[2] gives a minimum

plot(eval(alpha, data), h=-4..20, gridlines=true)

 
> 

# but all of this would have been much simpler if one had instanciated alpha
# before any calculus:

beta := eval(alpha, data);
fsolve(diff(beta, h));
eval(diff(beta, h$2), h=%);

beta := (2*((3/2)*piecewise(h < 1, 1, 2)+(piecewise(h < 1, 3, 4)+10)*(3+(3/2)*h/(5-h))))/(h*(piecewise(h < 1, 1, 2)+piecewise(h < 1, 3, 4)))

  

2.970655771

  

2.819541865

 (6)
> 

# note 1: this command is useful to capture discontinuities

discont(beta, h)

{0, 1, 5}

 (7)
> 

# minimize does not accept ranges of the form RealRange(Open(0), Open(5))
# and has no "avoid" option like in fsolve.

minimize(beta, h=0..5, location=true);

# Thus minimize cannot find the minimum unless you help it
minimize(beta, h=0.1..4.9, location=true);

-infinity, {[{h = 5}, -infinity]}

  

8.498780306, {[{h = 2.970655771}, 8.498780306]}

 (8)
> 

 


 

Download WithSolve.mw

Here is a possibility (one could write a...

mmcdara 7891
November 22 2020
0 3

Here is an (incomplete, see below) possibility.

The algorithm is given here  Greedy_algorithm_for_Egyptian_fractions

 

> 

restart:
> 

EgyptianFraction := proc(a)
  local b, c, s:
  b := 1/ceil(1/a);
  c := a-b:
  s := b:
  while c > 0 do
    b := 1/ceil(1/c);
    c := c-b:
    s := s, b:
  end do;
  return s
end proc:
> 

EgyptianFraction(7/15)

1/3, 1/8, 1/120

 (1)
> 

EgyptianFraction(5/121)

1/25, 1/757, 1/763309, 1/873960180913, 1/1527612795642093418846225

 (2)
> 

 


 

Download EgyptianFraction.mw


Watchout, the algorithm does not return necessarily what you had in mind.
For rationals > 1 the first fraction is 1/1, for instance

EgyptianFraction(13/12)
                                1 
                             1, --
                                12
# Expected 1/2 , 1/3 , 1/4 ?

To handle the case "rational > 1" thMAPLE code must be modified according to  Greedy_algorithm_for_Egyptian_fractions section "RELATED EXTENSIONS" (imposing the denominators to be > 1)

Another possibility is to use "brute force":

 

> 

restart:
> 

BruteForce := proc(a)
  local A, b, s:
  A := copy(a):

  s := NULL:
  b := 2;
  while A > 0 do
    while(A-1/b) < 0 do
      b := b+1:
    end do;
    A := A-1/b:
    s := s, 1/b:
    b := b+1:
  end do:
  return s;
end proc:
> 

BruteForce(13/12)

1/2, 1/3, 1/4

 (1)
> 

BruteForce(2)

1/2, 1/3, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9, 1/10, 1/15, 1/230, 1/57960

 (2)
> 

# do not try this one
# BruteForce(5/121)


 

Download BruteForce.mw

I tried to remain simple and pedagogic i...

mmcdara 7891
November 20 2020
0 1

I tried to remain simple and pedagogic in the attached file.
It tries to answer several points you raised in your post, even if I'm not sure I understood them correctly. 

As a rule some MAPLE's outputs are not reproduced on this site (for instance those which are generated by the Explore command); so do not pay too much attention to what is displayed above but load the attached file and execute into your own session.

Hope this will help you; and welcome on Mapleprimes

> 

restart:
> 

s := u*t + 1/2*a*t^2;

u*t+(1/2)*a*t^2

 (1)
> 

plots[interactive]();


What I understand by "I'd like to evaluate, graph, and table this"

> 

# 1/ evaluation for a given value of t (?)

eval(s, t=2);
#or
subs(t=2, s);
 

2*u+2*a

  

2*u+2*a

 (2)
> 

# 2/ graph
#    You have 3 symbolic quantities (a, t, u)
#
# Maybe you could begin by clicking on this Using the Exploration Assistant
#
# Modus operandi:
#   1/ On the main menu select Tools > Assistants > Plot Builder
#   2/ In the pop-up window select button <add>, type s in the pop-up window and <accept>
#      The assistant recognizes the expression of s and finds the three unknowns a, t, u
#      select <OK> to close this window
#   3/ In the new pop-up window:
#      3.1/ select the type of plot you want (upper combo box)
#           for instance <Interactive Plot with 2 parameters>
#      3.2/ choose the <x Axis> (for instance t) and define the plotting range
#      3.3/ for the 2 remaining parameters (here a and u) choose de ranges of variations
#      3.4/ select <Plot>
#
# On the plot you have 2 sliders you can move to "explore" the expression of s
#
# PS: running the assistant as I did automatically generates the command plots[interactive]();
# in the worksheet.
# Once faùmiliar with this assistant you will be capable to type this command and run it
# by yourself

 
> 

plots[interactive]();
> 

# A little bit more interesting
#
# Do not use the Assistant but type directly the line cmd := plots[interactive]();
#
# Instead of choosing <Plot> the last pop-up window, choose <Option>
# This will be help you to change a lot of plotting options.
# You can select <Plot> again, but select <Command> instead; you will get this

cmd := plots[interactive]();

Explore(plot(u*t+(1/2)*a*t^2, t = -10 .. 10, labels = [t, ""]), 'parameters' = [a = 0. .. 10., u = 0. .. 10.], 'initialvalues' = [a = 5.000000000, u = 5.000000000])

 (3)
> 

# You can now see exactly what MAPLE have understood and redo this same command

cmd
> 

# I guess that people familiar with MAPLE use directlys the explore command.
# Just look to the help page and look here Explore example worksheet  for several examples
# (you can copy-paste into your own worksheet the code chunks that suit you)
 
> 

# 3/ table ...
#    ... I'm not sure to understand correctly what you really want...
#    ... so tka this as a simple proposal
#
# First thing is that <table> is a particular structure in MAPLE.
# So, instead of using "tables", I will use a 2 columns matrix with
# the first containing different values of t and the second the corresponding
# values of s for given values of a and u.
#

# step 1: build a particular expression of s
#         Note this step is required only if you want a table of numerical values

specific_s := eval(s, [a=1, u=-1]);

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

 (4)
> 

# step 2: define the values of t, for instance between 0 and 3 by 0.5

t_values := [ seq(0..3, 0.5) ];  # this is just a possible way, among many, to do this

[0, .5, 1.0, 1.5, 2.0, 2.5, 3.0]

 (5)
> 

# step 3: evaluate <specific_s> for thes values of t
#         Here again there are a lot of ways to do this, for instance

s_values := [ seq( eval(specific_s, t=tau), tau in t_values) ]

[0, -.3750000000, -.5000000000, -.375000000, 0., .625000000, 1.500000000]

 (6)
> 

# step 3(bis): but what I prefer is to build an "arrow function" F : t ---> "the value of specific_s for this t"

F := unapply(specific_s, t)

proc (t) options operator, arrow; -t+(1/2)*t^2 end proc

 (7)
> 

# step 3(bis) continued: just do this (look in the help pages the "tilde" operator)

s_values := F~(t_values)

[0, -.3750000000, -.5000000000, -.375000000, 0., .625000000, 1.500000000]

 (8)
> 

# step 4: store t_values and s_values in a matrix

result := convert([t_values, s_values], Matrix);

# and transpose if you prefer

result := convert([t_values, s_values], Matrix)^+

result := Matrix(2, 7, {(1, 1) = 0, (1, 2) = .5, (1, 3) = 1.0, (1, 4) = 1.5, (1, 5) = 2.0, (1, 6) = 2.5, (1, 7) = 3.0, (2, 1) = 0, (2, 2) = -.3750000000, (2, 3) = -.5000000000, (2, 4) = -.375000000, (2, 5) = 0., (2, 6) = .625000000, (2, 7) = 1.500000000})

  

result := Matrix(7, 2, {(1, 1) = 0, (1, 2) = 0, (2, 1) = .5, (2, 2) = -.3750000000, (3, 1) = 1.0, (3, 2) = -.5000000000, (4, 1) = 1.5, (4, 2) = -.375000000, (5, 1) = 2.0, (5, 2) = 0., (6, 1) = 2.5, (6, 2) = .625000000, (7, 1) = 3.0, (7, 2) = 1.500000000})

 (9)
> 

# As I said there are many ways to do this, some are more efficient than others.
# If other Mapleprimes' members do reply to your question, they will surely give
# you different solutions.
#
# Once this matrix is build you can save it in a file (text file, binary file, ...
# and even xls/xlsx files).
# Just seek to the words <save>, <write>, <FileTools>, <Export>, <ExcelTools> and <Spread>
# in the help pages


About the maximum...

> 

# I understand that you already know of a particular value t=(v-u)/a of t?
#
# What you can do is evaluate F for this particular value of t

T := (v-u)/a;

F(T)

(v-u)/a

  

-(v-u)/a+(1/2)*(v-u)^2/a^2

 (10)
> 

# If you want to explore this new expression just run the command plots[interactive]();
# and proceed as previously.


 

Download A_little_help.mw

  I do not think MAPLE is capable ...

mmcdara 7891
November 20 2020
2 6

 

I do not think MAPLE is capable to solve such a complicated integral equation.

Neverteless you can always write a simple integration scheme to see what happens.
This is what I've done in the attached file. The work is anything but complete: the (elementary) numerical scheme has the lowest possible order (rectangle rule integration) ans some approximations have been done in order to simpligy the equation.

Nevertheless it could serve as a base for you to go further.
I hope someone else here will attack this problem of yours with more skill.

For a more correct resolution we need:

  1. the initial value of rho
  2. the maximum value of t you want to solve to
  3. an idea of some characteristic time in order to choose correctly the integration time step dt (different choices are presented)
     
> 

restart:
> 

w[c]:=0.01:delta:=5:Omega:=5:w[0]:=1:gamma_0:=1:k[b]:=1:T:=100:nu[n]:=2*Pi*n*k[b]*T:
> 

F:=2*gamma_0*w[c]^2*exp(-w[c]*tau)*sin(delta/Omega*(sin(Omega*t)-sin(Omega*(t-tau)))+w[0]*tau):
> 

k1:=2*k[b]*T*w[c]^2*sum((w[c]*exp(-w[c]*tau)-abs(nu[n])*exp(-abs(nu[n])*tau))/(w[c]^2-nu[n]^2),n=-infinity..infinity):
> 

G:=cos(delta/Omega*(sin(Omega*t)-sin(Omega*(t-tau)))+w[0]*tau)*k1:
> 

eq:=diff(rho(t),t)+Int(G*rho(tau),tau=0..t)=1/2*Int(G-F,tau=0..t):
 
> 

S := op(-1, k1);

opS := op(1, S);

opS0 := eval(opS, n=0);

sum((0.1000000000e-1*exp(-0.1000000000e-1*tau)-628.3185308*abs(n)*exp(-628.3185308*abs(n)*tau))/(-394784.1762*n^2+0.1000000000e-3), n = -infinity .. infinity)

  

(0.1000000000e-1*exp(-0.1000000000e-1*tau)-628.3185308*abs(n)*exp(-628.3185308*abs(n)*tau))/(-394784.1762*n^2+0.1000000000e-3)

  

100.0000000*exp(-0.1000000000e-1*tau)

 (1)
> 

# for n > 0 the denominator of opS is equivalent to -3.947841762*10^5*n^2 = (628.3185308*abs(n))^2 and the
# numerator is equivalent to 628.3185308*abs(n)*exp(-628.3185308*abs(n)*tau):
# Thus, assuming n > 0:

rel := Z = -628.3185308*n;

opS := Z*exp(Z*tau)/(-Z^2):
opS := eval(opS, rel);

Z = -628.3185308*n

  

0.1591549431e-2*exp(-628.3185308*n*tau)/n

 (2)
> 

S :=opS0 + sum(opS, n=1..+infinity)

100.0000000*exp(-0.1000000000e-1*tau)-0.1591549431e-2*ln(1.-1./exp(628.3185308*tau))

 (3)
> 

k1 := 2*k[b]*T*w[c]^2*S

2.000000000*exp(-0.1000000000e-1*tau)-0.3183098862e-4*ln(1.-1./exp(628.3185308*tau))

 (4)
> 

G:=cos(delta/Omega*(sin(Omega*t)-sin(Omega*(t-tau)))+w[0]*tau)*k1;

cos(sin(5*t)-sin(5*t-5*tau)+tau)*(2.000000000*exp(-0.1000000000e-1*tau)-0.3183098862e-4*ln(1.-1./exp(628.3185308*tau)))

 (5)
> 

# for tau >> 0 G is equivalent to

GL := (cos(sin(5*t)-sin(5*t-5*tau)+tau)*(2.000000000*exp(-0.1000000000e-1*tau)));

2.000000000*cos(sin(5*t)-sin(5*t-5*tau)+tau)*exp(-0.1000000000e-1*tau)

 (6)
> 

RHS := eval(-G*rho(tau)+1/2(G-F), rho(tau)=U);

-cos(sin(5*t)-sin(5*t-5*tau)+tau)*(2.000000000*exp(-0.1000000000e-1*tau)-0.3183098862e-4*ln(1.-1./exp(628.3185308*tau)))*U+1/2

 (7)
> 

dt    := 0.002:
tend  := 1:
times := [seq(0..tend, dt)]:
nstep := numelems(times):
sol   := Vector(nstep):

sol[1] := 0: # assuming an initial condition rho(0)=0

for n from 2 to nstep do
  sol[n] := sol[n-1] + dt * add( eval(RHS, [t=n*dt, tau=(m+1/2)*dt, U=sol[m]]), m=1..n )
end do:

plot([seq([times[n], sol[n]], n=1..nstep)])

 
> 

dt    := 0.001:
tend  := 1:
times := [seq(0..tend, dt)]:
nstep := numelems(times):
sol   := Vector(nstep):

sol[1] := 0: # assuming an initial condition rho(0)=0

for n from 2 to nstep do
  sol[n] := sol[n-1] + dt * add( eval(RHS, [t=n*dt, tau=(m+1/2)*dt, U=sol[m]]), m=1..n )
end do:

plot([seq([times[n], sol[n]], n=1..nstep)])

 
> 

dt    := 0.0005:
tend  := 1:
times := [seq(0..tend, dt)]:
nstep := numelems(times):
sol   := Vector(nstep):

sol[1] := 0: # assuming an initial condition rho(0)=0

for n from 2 to nstep do
  sol[n] := sol[n-1] + dt * add( eval(RHS, [t=n*dt, tau=(m+1/2)*dt, U=sol[m]]), m=1..n )
end do:

plot([seq([times[n], sol[n]], n=1..nstep)])

 
> 

 


 

Download abcde.mw

 

 

 

Hi, Here is a generic procedure to ...

mmcdara 7891
November 20 2020
3 5

Hi, 

Here is a generic procedure to build a system of equations like those you seek.
It is not restricted to three rooms (+1 for the exterior), nor to a single source term (=furnace), and some walls can be missing ("blind room" without communication to the exterior, adiabatic wall, rooms along a corridor, ...).

To be generic the tranmission coefficients (k1, k2, ...) are renamed k_{r, r'} where r and r' are two rooms; if no wall exist between r and r' then k_{r, r'} is set to 0.

BONUS: if you want to visualize the cubicle you can do this

with(GraphTheory):
cubicle := Graph({walls[]}):
DrawGraph(cubicle, style=`if`(IsPlanar(cubicle), planar, NULL))


 

> 

restart:
> 

interface(version)

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

 (1)
> 

GenSys := proc(rooms, walls, sources)
  local balance, nr, temp, direct_fluxes, reverse_fluxes, all_fluxes, all_exchanges:
  local null_k, equations, ics:

  balance := proc(i)
    local r, from_r_to_others, local_fluxes:
    r := rooms[i]:
    from_r_to_others := map(u -> if u[1]=r then u end if, all_exchanges):
    local_fluxes     := map(u -> eval(phi_[u], all_fluxes), from_r_to_others);

    # remove terms that contain fluxes through non existing walls

    null_k           := map(u -> if `not`(member(op(1,u), walls)) then u=0 end if, op~(1, local_fluxes)):
    local_fluxes     := eval(local_fluxes, null_k):

    if r in op~(1, sources) then
      diff(temp[i], t) = - add(local_fluxes) + select(has, sources, r)[][2]
    else
      diff(temp[i], t) = - add(local_fluxes)
    end if:
  end proc:

  nr     := numelems(rooms):
  temp   := seq(theta[rooms[i]](t), i=1..nr);

  direct_fluxes :=
        [
          seq(
            seq(
              phi_[[rooms[i], rooms[j]]] = k[{rooms[i], rooms[j]}]*(temp[i]-temp[j]),
              j=i+1..nr
            ),
            i=1..nr-1
          )
        ]:

  reverse_fluxes := map(u -> eval(u, op(lhs(u)) =~ ListTools:-Reverse(op(lhs(u)))), direct_fluxes):

  all_fluxes    := [direct_fluxes[], reverse_fluxes[]];
  all_exchanges := map(u -> op(1, lhs(u)), all_fluxes);

  equations := { seq(balance(i), i=1..nr-1) }:
  ics       := { seq(eval(temp[i], t=0) = Theta[rooms[i]], i=1..nr-1) }:

  return equations union ics
end proc:


STAGE 1: GENERATE THE DIFFERENTIAL SYSTEM

 

 

> 

# Ext is the "exterior room" and must be declared in the last position
#
# rooms   : list of the names of the rooms
# walls   : list of separating walls
#           each wall is a set {r, r'} where r and r' are two different rooms
# sources : list of source terms
#           each source term is a list [r, s] where r is a room name and s the expression of the source term
#
# Temperatures are denote by theta and initial temperatures by Theta
# To be general, the conduction coefficient through the wall that separates roons r and r' is
# is denoted k_{r, r'} (the braces mean that k_{r, r'} = k_{r', r}

rooms   := [A, B, C, Ext];
nr      := numelems(rooms):
walls   := [ seq(seq({rooms[i], rooms[j]}, j=i+1..nr), i=1..nr-1) ];
sources := [ [C, F(t)] ];

print();

sys   := GenSys(rooms, walls, sources):
print~(sys):

[A, B, C, Ext]

  

[{A, B}, {A, C}, {A, Ext}, {B, C}, {B, Ext}, {C, Ext}]

  

[[C, F(t)]]

  

 

diff(theta[A](t), t) = -k[{A, B}]*(theta[A](t)-theta[B](t))-k[{A, C}]*(theta[A](t)-theta[C](t))-k[{A, Ext}]*(theta[A](t)-theta[Ext](t))

  

diff(theta[B](t), t) = -k[{B, C}]*(theta[B](t)-theta[C](t))-k[{B, Ext}]*(theta[B](t)-theta[Ext](t))-k[{A, B}]*(theta[B](t)-theta[A](t))

  

diff(theta[C](t), t) = -k[{C, Ext}]*(theta[C](t)-theta[Ext](t))-k[{A, C}]*(theta[C](t)-theta[A](t))-k[{B, C}]*(theta[C](t)-theta[B](t))+F(t)

  

theta[A](0) = Theta[A]

  

theta[B](0) = Theta[B]

  

theta[C](0) = Theta[C]

 (2)
> 

# Exemple of a little bit more complicated example
# (more rooms, some walls missing, more heat sources)

if false then
  rooms   := [A, B, C, E, Ext]:
  nr      := numelems(rooms):
  walls   := [ seq(seq({rooms[i], rooms[j]}, j=min(nr, i+2)..nr), i=1..nr-1) ]:
  print(walls):
  sources := [ [C, F(t)], [A, g(t)] ]:

  sys   := GenSys(rooms, walls, sources):
  print~(sys):
end if:


STAGE 2: STATIONNARY STATE

 

> 

Unknowns := convert(op~(select(has, lhs~(sys), diff)) minus {t}, list)

[theta[A](t), theta[B](t), theta[C](t)]

 (3)
> 

# Equilibrium temperatures
#
# I guess "Equilibrium temperatures" means "stationnary solution"?
#
# L__r represents the stationnary temperature in room r

Stationnary_Unknowns := [seq(L__||r, r in rooms[1..-2])];


Stationnary_System := rhs~(sys[1..numelems(rooms)-1]):
Stationnary_System := eval(Stationnary_System, Unknowns =~ Stationnary_Unknowns):
print~(Stationnary_System):
 

[L__A, L__B, L__C]

  

-k[{A, B}]*(L__A-L__B)-k[{A, C}]*(L__A-L__C)-k[{A, Ext}]*(L__A-theta[Ext](t))

  

-k[{B, C}]*(L__B-L__C)-k[{B, Ext}]*(L__B-theta[Ext](t))-k[{A, B}]*(L__B-L__A)

  

-k[{C, Ext}]*(L__C-theta[Ext](t))-k[{A, C}]*(L__C-L__A)-k[{B, C}]*(L__C-L__B)+F(t)

 (4)
> 

# Simplified case
#   F(t) = S =constant
#   theta__Ext(t) = Theta__Ext =constant

Example := eval(Stationnary_System, {F(t)=S, theta[Ext](t)=Theta[Ext]}):
print~(Example):

-k[{A, B}]*(L__A-L__B)-k[{A, C}]*(L__A-L__C)-k[{A, Ext}]*(L__A-Theta[Ext])

  

-k[{B, C}]*(L__B-L__C)-k[{B, Ext}]*(L__B-Theta[Ext])-k[{A, B}]*(L__B-L__A)

  

-k[{C, Ext}]*(L__C-Theta[Ext])-k[{A, C}]*(L__C-L__A)-k[{B, C}]*(L__C-L__B)+S

 (5)
> 

# MAPLE can solve this system but the solution is very lengthy

Stationnary_Solution := solve(Example, Stationnary_Unknowns):
length(%);

6762

 (6)
> 

MyRules := [
             k[{A, B}]   = k__3,
             k[{A, C}]   = k__3,
             k[{A, Ext}] = k__4,
             k[{B, C}]   = k__1,
             k[{B, Ext}] = k__2,
             k[{C, Ext}] = k__5
           ]:

collect(simplify(eval(L__A, Stationnary_Solution[1][1])), [Theta[Ext], k[{A, B}], k_[{A, C}], k_[{A, Ext}]]):
eval(%, MyRules)

 

Theta[Ext]+((S*k__1+S*k__3)*k__3+S*k__3*k__1+S*k__3*k__2)/((k__1*k__2+k__1*k__4+k__1*k__5+k__2*k__3+k__2*k__5+k__3*k__4+k__3*k__5+k__4*k__5)*k__3+k__3*k__4*k__1+k__3*k__4*k__2+k__3*k__1*k__2+k__3*k__1*k__5+k__3*k__2*k__5+k__4*k__1*k__2+k__4*k__1*k__5+k__4*k__2*k__5)

 (7)
> 

 

 


 

Download General_Model.mw

Hi, I think your text contains unne...

mmcdara 7891
November 19 2020
1 1

Hi, 

I think your text contains unnecessary material and that your question simply resume to "How can I test if a sample can be considered as a sample of a given random variable?"

This is a recurring question on this site and has been answered regularly.
Here the originality comes from the fact the WignerGOE distribution you refer to is not implemented in Maple.
This issue can be easily by-passed as described in the attached file.


 

> 

restart:
> 

with(Statistics):
> 

# f is assumed to be the pdf of some random variable W (=WignerGOE)

f := n -> Pi/2*n*exp(-Pi/4*n^2)

proc (n) options operator, arrow; (1/2)*Pi*n*exp(-(1/4)*Pi*n^2) end proc

 (1)
> 

# let's check this

int(f(n), n=+infinity..+infinity);  # obvious for f(n) is odd

0

 (2)
> 

# obviously f(n) is not a pdf ...
# ... butf(n) restricted to n>=0 is one:

int(f(n), n=0..+infinity)

1

 (3)
> 

# redefine f(n) correctly:

pdf := unapply(piecewise(n<0, 0, f(n)), n);

proc (n) options operator, arrow; piecewise(n < 0, 0, (1/2)*Pi*n*exp(-(1/4)*Pi*n^2)) end proc

 (4)
> 

# for future needs it could be useful to define also a cdf:

cdf := unapply(int(f(z), z=0..n), n);

proc (n) options operator, arrow; 1-exp(-(1/4)*Pi*n^2) end proc

 (5)
> 

# use pdf(n) to define a new random variable

W := Distribution(PDF = pdf, CDF = cdf)

_m4724647648

 (6)
> 

# example of use

CodeTools:-Usage( Histogram(Sample(W, 10^3, method='envelope')) )

memory used=60.17MiB, alloc change=46.01MiB, cpu time=759.00ms, real time=760.00ms, gc time=37.25ms

  

 
> 

# Generate now a sample from another distribution which looks like
# the one above.
# One may think to a Generalized Beta distribution of the form c * Beta(a, b)
#
# First step : compute the mean and variance of W

mv__W := [ (Mean, Variance)(W) ]

[(1/2)*4^(1/2), -(Pi-4)/Pi]

 (7)
> 

# second step : determine the values of the parameters a and b such that the
#               mean and variance of the Generalized Beta equal those of W

GenBeta := c*RandomVariable(BetaDistribution(a, b)):
mv__GB  := [ (Mean, Variance)(GenBeta) ];
sol     := solve(mv__W =~ mv__GB, [a, b]);

# Choose a positive value for c: it will be large enough for the support of
# the generalized Beta distribution encompasses the W's (wich is infinite).
# For instance:

C     := 5:
SOL   := evalf(eval(sol, c=C))[];

FittedGenBeta := C*RandomVariable(BetaDistribution(op(rhs~(SOL))));

[c*a/(a+b), c^2*a*b/((a+b)^2*(a+b+1))]

  

[[a = -(Pi*c-4)/(c*(Pi-4)), b = -(Pi*c^2-Pi*c-4*c+4)/(c*(Pi-4))]]

  

[a = 2.727833894, b = 10.91133558]

  

5*_R3

 (8)
> 

# Now draw a sample from FittedGenBeta and compare its histogram to the pdf of W:

S__FGB := Sample(FittedGenBeta, 10^4):

infolevel[Statistics] := 0:
plots:-display(
  Histogram(S__FGB, minbins=30, range=0..C, view=[0..C, default], transparency=0.3),
  plot(pdf(n), n=0..C, color=red,thickness=3)
)

 
> 

# Now this is, not the Wigner's surmise but MY hypothesis:
#
# H0 : "the sample S is drawn from the random variable W"
#
# Let's check it

infolevel[Statistics] := 1:
ChiSquareSuitableModelTest(S__FGB, W, bins=20):
 

Chi-Square Test for Suitable Probability Model
----------------------------------------------
Null Hypothesis:
Sample was drawn from specified probability distribution
Alt. Hypothesis:
Sample was not drawn from specified probability distribution
Bins:                    20
Degrees of freedom:      19
Distribution:            ChiSquare(19)
Computed statistic:      55.004
Computed pvalue:         2.3211e-05
Critical value:          30.1435273589987
Result: [Rejected]
This statistical test provides evidence that the null hypothesis is false

  
> 

# To be sure that ChiSquareSuitableModelTest works correctly,
# let's verify what this this gives when the sample is really drawn from W:

S__W := Sample(W, 10^4, method='envelope'):
ChiSquareSuitableModelTest(S__W, W, bins=20):

Chi-Square Test for Suitable Probability Model
----------------------------------------------
Null Hypothesis:
Sample was drawn from specified probability distribution
Alt. Hypothesis:
Sample was not drawn from specified probability distribution
Bins:                    20
Degrees of freedom:      19
Distribution:            ChiSquare(19)
Computed statistic:      14.776
Computed pvalue:         0.736725
Critical value:          30.1435273589987
Result: [Accepted]
This statistical test does not provide enough evidence to conclude that the null hypothesis is false

  
> 

 


 

Download WignerGOE_2.mw

Here is a way which uses the full power ...

mmcdara 7891
November 17 2020
1 0

Here is a way which uses the full power of Maple in formal manipulation of random variables
This in my oxn opinion the best way to proceed.

restart:
with(Statistics):

# probabilies for each event
P := [$1..4] /~ 10:

# define 4 iid random variables whose probabilities are given by P
for m from 1 to 4 do
  X||m := RandomVariable(ProbabilityTable(P))
end do:

# define new RV as the sum of several (from 2 to 4) "basic" RVs
for n from 2 to 4 do
  S||n := add(X||m, m=1..n)
end do:

# a few examples

for s from 2 to 8 do
  printf("The probability that the sum of 2 RVs equals %d is %a\n", s, Probability(S2=s))
end do:

for s from 3 to 12 do
  printf("The probability that the sum of 3 RVs equals %d is %a\n", s, Probability(S3=s))
end do:

Download A_Way.mw

If you want something more "computational" you can do this (here is a particular example
 

# prob to get 10 with the sum of 3 RVs

[ map(i -> if max(i) < 5 then i end if, combinat:-composition(10, 3))[] ];
map(i -> i/~10, %);
map(mul, %);
add(%)
51 
---
250

 

There exist a lot of methods to compare ...

mmcdara 7891
November 15 2020
0 0

There exist a lot of methods to compare an empirical distribution (by abusive language an histogram)  with a theoritical distribution.

One of the most classical method is the Chi-Square Goodness of Fit (Chi 2 GOF for short)
It's accessible in MAPLE via  the Statistics package, look at the help page Statistics[ChiSquareSuitableModelTest].
Nevertheless the way this test is implemented is questionnable (I've akready post something about that).

I also remember that @Carl Love  published a post about the G-test (wich can be saw as a variant of the Chi GOF test).

In the particular case where you want to test if your empirical distribution can be considered as a gaussian, you enter the class of normality tests.
Here again there exist plenty of them, for instance the Shapiro-Wilk test (see Statistics[ShapiroWilkWTest]).

Let me know if you meet any difficulties to use these tests.

Here are a few ideas: use "paramete...

mmcdara 7891
November 09 2020
0 0

Here are a few ideas:

  • use "parameters" to say MAPLE that xlow is a parameter you will fix later,
  • use try..catch..end try for your ode presents a singularity at some time depending on the value of slow you use
    (lasterror is used to capture the time at which this singularity appears (if any)),
  • ue plottools:-transform to exchange the role of x and y
     

Hope you will find this useful

PS: a prori there is no need to use the option range unless tou want to refine the plot.

> 

restart
> 

ode := diff(x(t), t) = 16250.25391*(1 - (487*x(t))/168 + 4*Pi*x(t)^(3/2) + (274229*x(t)^2)/72576 - (254*Pi*x(t)^(5/2))/21 + (119.6109573 - (856*ln(16*x(t)))/105)*x(t)^3 + (30295*Pi*x(t)^(7/2))/1728 + (7.617741607 - 23.53000000*ln(x(t)))*x(t)^4 + (535.2001594 - 102.446*ln(x(t)))*x(t)^(9/2) + (413.8828821 + 323.5521650*ln(x(t)))*x(t)^5 + (1533.899179 - 390.2690000*ln(x(t)))*x(t)^(11/2) + (2082.250556 + 423.6762500*ln(x(t)) + 33.2307*ln(x(t)^2))*x(t)^6)*x(t)^5;
> 

ics  := x(0) = xlow;
tfin := 10;

solx := dsolve({ics, ode}, numeric, parameters=[xlow], method=rosenbrock);
solx(parameters=[xlow=1e-1]):

try
  solx(tfin):
  p := plots[odeplot](solx, 0 .. tmax):
catch:
  print(lasterror);
  tmax := lasterror()[-1];
  p := plots[odeplot](solx, 0 .. tmax):
end try:

plots[display](p)

diff(x(t), t) = 16250.25391*(1-(487/168)*x(t)+4*Pi*x(t)^(3/2)+(274229/72576)*x(t)^2-(254/21)*Pi*x(t)^(5/2)+(119.6109573-(856/105)*ln(16*x(t)))*x(t)^3+(30295/1728)*Pi*x(t)^(7/2)+(7.617741607-23.53000000*ln(x(t)))*x(t)^4+(535.2001594-102.446*ln(x(t)))*x(t)^(9/2)+(413.8828821+323.5521650*ln(x(t)))*x(t)^5+(1533.899179-390.2690000*ln(x(t)))*x(t)^(11/2)+(2082.250556+423.6762500*ln(x(t))+33.2307*ln(x(t)^2))*x(t)^6)*x(t)^5

  

x(0) = xlow

  

10

  

proc (x_rosenbrock) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if 1 < nargs then error "invalid input: too many arguments" end if; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then _xout := evalf[_EnvDSNumericSaveDigits](x_rosenbrock) else _xout := evalf(x_rosenbrock) end if; _dat := Array(1..4, {(1) = proc (_xin) local _xout, _dtbl, _dat, _vmap, _x0, _y0, _val, _dig, _n, _ne, _nd, _nv, _pars, _ini, _par, _i, _j, _k, _src; option `Copyright (c) 2002 by Waterloo Maple Inc. All rights reserved.`; table( [( "complex" ) = false ] ) _xout := _xin; _pars := [xlow = xlow]; _dtbl := array( 1 .. 4, [( 1 ) = (array( 1 .. 24, [( 1 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 2 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 3 ) = ([0, 0, 0, Array(1..0, 1..2, {}, datatype = float[8], order = C_order)]), ( 4 ) = (Array(1..54, {(1) = 1, (2) = 1, (3) = 0, (4) = 0, (5) = 1, (6) = 0, (7) = 0, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 0, (14) = 0, (15) = 0, (16) = 0, (17) = 0, (18) = 0, (19) = 30000, (20) = 0, (21) = 0, (22) = 2, (23) = 3, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 2, (34) = 0, (35) = 0, (36) = 0, (37) = 0, (38) = 0, (39) = 0, (40) = 0, (41) = 0, (42) = 0, (43) = 1, (44) = 0, (45) = 0, (46) = 0, (47) = 0, (48) = 0, (49) = 0, (50) = 50, (51) = 1, (52) = 0, (53) = 0, (54) = 0}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = .0, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = .0, (6) = .0, (7) = .0, (8) = 0.10e-5, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = 1.0, (14) = .0, (15) = .49999999999999, (16) = .0, (17) = 1.0, (18) = 1.0, (19) = .0, (20) = .0, (21) = 1.0, (22) = 1.0, (23) = .0, (24) = .0, (25) = 0.10e-14, (26) = .0, (27) = .0, (28) = .0}, datatype = float[8], order = C_order)), ( 6 ) = (Array(1..2, {(1) = xlow, (2) = Float(undefined)})), ( 7 ) = ([Array(1..4, 1..7, {(1, 1) = .0, (1, 2) = .203125, (1, 3) = .3046875, (1, 4) = .75, (1, 5) = .8125, (1, 6) = .40625, (1, 7) = .8125, (2, 1) = 0.6378173828125e-1, (2, 2) = .0, (2, 3) = .279296875, (2, 4) = .27237892150878906, (2, 5) = -0.9686851501464844e-1, (2, 6) = 0.1956939697265625e-1, (2, 7) = .5381584167480469, (3, 1) = 0.31890869140625e-1, (3, 2) = .0, (3, 3) = -.34375, (3, 4) = -.335235595703125, (3, 5) = .2296142578125, (3, 6) = .41748046875, (3, 7) = 11.480712890625, (4, 1) = 0.9710520505905151e-1, (4, 2) = .0, (4, 3) = .40350341796875, (4, 4) = 0.20297467708587646e-1, (4, 5) = -0.6054282188415527e-2, (4, 6) = -0.4770040512084961e-1, (4, 7) = .77858567237854}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = 1.0, (2, 1) = .25, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = 1.0, (3, 1) = .1875, (3, 2) = .5625, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = 2.0, (4, 1) = .23583984375, (4, 2) = -.87890625, (4, 3) = .890625, (4, 4) = .0, (4, 5) = .0, (4, 6) = .2681884765625, (5, 1) = .1272735595703125, (5, 2) = -.5009765625, (5, 3) = .44921875, (5, 4) = -0.128936767578125e-1, (5, 5) = .0, (5, 6) = 0.626220703125e-1, (6, 1) = -0.927734375e-1, (6, 2) = .626220703125, (6, 3) = -.4326171875, (6, 4) = .1418304443359375, (6, 5) = -0.861053466796875e-1, (6, 6) = .3131103515625}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .386, (3) = .21, (4) = .63, (5) = 1.0, (6) = 1.0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .25, (2) = -.1043, (3) = .1035, (4) = -0.362e-1, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 1.544, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = .9466785280815533, (3, 2) = .25570116989825814, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = 3.3148251870684886, (4, 2) = 2.896124015972123, (4, 3) = .9986419139977808, (4, 4) = .0, (4, 5) = .0, (5, 1) = 1.2212245092262748, (5, 2) = 6.019134481287752, (5, 3) = 12.537083329320874, (5, 4) = -.687886036105895, (5, 5) = .0, (6, 1) = 1.2212245092262748, (6, 2) = 6.019134481287752, (6, 3) = 12.537083329320874, (6, 4) = -.687886036105895, (6, 5) = 1.0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = -5.6688, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = -2.4300933568337584, (3, 2) = -.20635991570891224, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = -.10735290581452621, (4, 2) = -9.594562251021896, (4, 3) = -20.470286148096154, (4, 4) = .0, (4, 5) = .0, (5, 1) = 7.496443313968615, (5, 2) = -10.246804314641219, (5, 3) = -33.99990352819906, (5, 4) = 11.708908932061595, (5, 5) = .0, (6, 1) = 8.083246795922411, (6, 2) = -7.981132988062785, (6, 3) = -31.52159432874373, (6, 4) = 16.319305431231363, (6, 5) = -6.0588182388340535}, datatype = float[8], order = C_order), Array(1..3, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 10.126235083446911, (2, 2) = -7.487995877607633, (2, 3) = -34.800918615557414, (2, 4) = -7.9927717075687275, (2, 5) = 1.0251377232956207, (3, 1) = -.6762803392806898, (3, 2) = 6.087714651678606, (3, 3) = 16.43084320892463, (3, 4) = 24.767225114183653, (3, 5) = -6.5943891257167815}, datatype = float[8], order = C_order)]), ( 9 ) = ([Array(1..1, {(1) = 1.0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = 0}, datatype = integer[8]), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order)]), ( 8 ) = ([Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..1, {(1, 1) = .0, (2, 0) = .0, (2, 1) = .0, (3, 0) = .0, (3, 1) = .0, (4, 0) = .0, (4, 1) = .0, (5, 0) = .0, (5, 1) = .0, (6, 0) = .0, (6, 1) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = x(t)]`; if Y[1] < 0 then YP[1] := undefined; return 0 end if; if Y[1] < 0 then YP[1] := undefined; return 0 end if; if Y[1] < 0 then YP[1] := undefined; return 0 end if; if Y[1] < 0 then YP[1] := undefined; return 0 end if; if Y[1] < 0 then YP[1] := undefined; return 0 end if; YP[1] := 16250.25391*(1-(487/168)*Y[1]+12.5663706143592*evalf(Y[1]^(3/2))+(274229/72576)*Y[1]^2-37.9983111434194*evalf(Y[1]^(5/2))+(119.6109573-(856/105)*ln(16*Y[1]))*Y[1]^3+55.0778642595502*evalf(Y[1]^(7/2))+(7.617741607-23.53000000*ln(Y[1]))*Y[1]^4+(535.2001594-102.446*ln(Y[1]))*evalf(Y[1]^(9/2))+(413.8828821+323.5521650*ln(Y[1]))*Y[1]^5+(1533.899179-390.2690000*ln(Y[1]))*evalf(Y[1]^(11/2))+(2082.250556+423.6762500*ln(Y[1])+33.2307*ln(Y[1]^2))*Y[1]^6)*Y[1]^5; 0 end proc, proc (X, Y, FX, FY) FX[1 .. 1] := 0; FY[1 .. 1, 1 .. 1] := 0; if Y[1] < 0 then FX[1] := undefined; return 0 end if; if Y[1] < 0 then FX[1] := undefined; return 0 end if; if Y[1] < 0 then FX[1] := undefined; return 0 end if; if Y[1] < 0 then FX[1] := undefined; return 0 end if; if Y[1] < 0 then FX[1] := undefined; return 0 end if; if Y[1] < 0 then FX[1] := undefined; return 0 end if; FY[1, 1] := 16250.25391*(-487/168+18.8495559215388*evalf(Y[1]^(1/2))+(274229/36288)*Y[1]-94.9957778585485*evalf(Y[1]^(3/2))-(856/105)*Y[1]^2+3*(119.6109573-(856/105)*ln(16*Y[1]))*Y[1]^2+192.772524908426*evalf(Y[1]^(5/2))-23.53000000*Y[1]^3+4*(7.617741607-23.53000000*ln(Y[1]))*Y[1]^3-102.446*evalf(Y[1]^(7/2))+(9/2)*(535.2001594-102.446*ln(Y[1]))*evalf(Y[1]^(7/2))+323.5521650*Y[1]^4+5*(413.8828821+323.5521650*ln(Y[1]))*Y[1]^4-390.2690000*evalf(Y[1]^(9/2))+(11/2)*(1533.899179-390.2690000*ln(Y[1]))*evalf(Y[1]^(9/2))+490.1376500*Y[1]^5+6*(2082.250556+423.6762500*ln(Y[1])+33.2307*ln(Y[1]^2))*Y[1]^5)*Y[1]^5+81251.26955*(1-(487/168)*Y[1]+12.5663706143592*evalf(Y[1]^(3/2))+(274229/72576)*Y[1]^2-37.9983111434194*evalf(Y[1]^(5/2))+(119.6109573-(856/105)*ln(16*Y[1]))*Y[1]^3+55.0778642595502*evalf(Y[1]^(7/2))+(7.617741607-23.53000000*ln(Y[1]))*Y[1]^4+(535.2001594-102.446*ln(Y[1]))*evalf(Y[1]^(9/2))+(413.8828821+323.5521650*ln(Y[1]))*Y[1]^5+(1533.899179-390.2690000*ln(Y[1]))*evalf(Y[1]^(11/2))+(2082.250556+423.6762500*ln(Y[1])+33.2307*ln(Y[1]^2))*Y[1]^6)*Y[1]^4; 0 end proc, 0, 0, 0, 0, 0, 0]), ( 13 ) = (), ( 12 ) = (), ( 15 ) = ("rosenbrock"), ( 14 ) = ([0, 0]), ( 18 ) = ([]), ( 19 ) = (0), ( 16 ) = ([0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) option `[Y[1] = x(t)]`; if Y[1] < 0 then YP[1] := undefined; return 0 end if; if Y[1] < 0 then YP[1] := undefined; return 0 end if; if Y[1] < 0 then YP[1] := undefined; return 0 end if; if Y[1] < 0 then YP[1] := undefined; return 0 end if; if Y[1] < 0 then YP[1] := undefined; return 0 end if; YP[1] := 16250.25391*(1-(487/168)*Y[1]+12.5663706143592*evalf(Y[1]^(3/2))+(274229/72576)*Y[1]^2-37.9983111434194*evalf(Y[1]^(5/2))+(119.6109573-(856/105)*ln(16*Y[1]))*Y[1]^3+55.0778642595502*evalf(Y[1]^(7/2))+(7.617741607-23.53000000*ln(Y[1]))*Y[1]^4+(535.2001594-102.446*ln(Y[1]))*evalf(Y[1]^(9/2))+(413.8828821+323.5521650*ln(Y[1]))*Y[1]^5+(1533.899179-390.2690000*ln(Y[1]))*evalf(Y[1]^(11/2))+(2082.250556+423.6762500*ln(Y[1])+33.2307*ln(Y[1]^2))*Y[1]^6)*Y[1]^5; 0 end proc, proc (X, Y, FX, FY) FX[1 .. 1] := 0; FY[1 .. 1, 1 .. 1] := 0; if Y[1] < 0 then FX[1] := undefined; return 0 end if; if Y[1] < 0 then FX[1] := undefined; return 0 end if; if Y[1] < 0 then FX[1] := undefined; return 0 end if; if Y[1] < 0 then FX[1] := undefined; return 0 end if; if Y[1] < 0 then FX[1] := undefined; return 0 end if; if Y[1] < 0 then FX[1] := undefined; return 0 end if; FY[1, 1] := 16250.25391*(-487/168+18.8495559215388*evalf(Y[1]^(1/2))+(274229/36288)*Y[1]-94.9957778585485*evalf(Y[1]^(3/2))-(856/105)*Y[1]^2+3*(119.6109573-(856/105)*ln(16*Y[1]))*Y[1]^2+192.772524908426*evalf(Y[1]^(5/2))-23.53000000*Y[1]^3+4*(7.617741607-23.53000000*ln(Y[1]))*Y[1]^3-102.446*evalf(Y[1]^(7/2))+(9/2)*(535.2001594-102.446*ln(Y[1]))*evalf(Y[1]^(7/2))+323.5521650*Y[1]^4+5*(413.8828821+323.5521650*ln(Y[1]))*Y[1]^4-390.2690000*evalf(Y[1]^(9/2))+(11/2)*(1533.899179-390.2690000*ln(Y[1]))*evalf(Y[1]^(9/2))+490.1376500*Y[1]^5+6*(2082.250556+423.6762500*ln(Y[1])+33.2307*ln(Y[1]^2))*Y[1]^5)*Y[1]^5+81251.26955*(1-(487/168)*Y[1]+12.5663706143592*evalf(Y[1]^(3/2))+(274229/72576)*Y[1]^2-37.9983111434194*evalf(Y[1]^(5/2))+(119.6109573-(856/105)*ln(16*Y[1]))*Y[1]^3+55.0778642595502*evalf(Y[1]^(7/2))+(7.617741607-23.53000000*ln(Y[1]))*Y[1]^4+(535.2001594-102.446*ln(Y[1]))*evalf(Y[1]^(9/2))+(413.8828821+323.5521650*ln(Y[1]))*Y[1]^5+(1533.899179-390.2690000*ln(Y[1]))*evalf(Y[1]^(11/2))+(2082.250556+423.6762500*ln(Y[1])+33.2307*ln(Y[1]^2))*Y[1]^6)*Y[1]^4; 0 end proc, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 24 ) = (0) ] )) ] ); _y0 := Array(0..2, {(1) = 0., (2) = xlow}); _vmap := array( 1 .. 1, [( 1 ) = (1) ] ); _x0 := _dtbl[1][5][5]; _n := _dtbl[1][4][1]; _ne := _dtbl[1][4][3]; _nd := _dtbl[1][4][4]; _nv := _dtbl[1][4][16]; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then if _Env_smart_dsolve_numeric = true or _dtbl[1][4][10] = 1 then if _xout = "left" then if type(_dtbl[2], 'table') then return _dtbl[2][5][1] end if elif _xout = "right" then if type(_dtbl[3], 'table') then return _dtbl[3][5][1] end if end if end if; return _dtbl[1][5][5] elif _xout = "method" then return _dtbl[1][15] elif _xout = "storage" then return evalb(_dtbl[1][4][10] = 1) elif _xout = "leftdata" then if not type(_dtbl[2], 'array') then return NULL else return eval(_dtbl[2]) end if elif _xout = "rightdata" then if not type(_dtbl[3], 'array') then return NULL else return eval(_dtbl[3]) end if elif _xout = "enginedata" then return eval(_dtbl[1]) elif _xout = "enginereset" then _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); return NULL elif _xout = "initial" then return procname(_y0[0]) elif _xout = "laxtol" then return _dtbl[`if`(member(_dtbl[4], {2, 3}), _dtbl[4], 1)][5][18] elif _xout = "numfun" then return `if`(member(_dtbl[4], {2, 3}), _dtbl[_dtbl[4]][4][18], 0) elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return procname(_y0[0]), [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _dtbl[4] <> 2 and _dtbl[4] <> 3 or _x0-_dtbl[_dtbl[4]][5][1] = 0. then error "no information is available on last computed point" else _xout := _dtbl[_dtbl[4]][5][1] end if elif _xout = "function" then if _dtbl[1][4][33]-2. = 0 then return eval(_dtbl[1][10], 1) else return eval(_dtbl[1][10][1], 1) end if elif _xout = "map" then return copy(_vmap) elif type(_xin, `=`) and type(rhs(_xin), 'list') and member(lhs(_xin), {"initial", "parameters", "initial_and_parameters"}) then _ini, _par := [], []; if lhs(_xin) = "initial" then _ini := rhs(_xin) elif lhs(_xin) = "parameters" then _par := rhs(_xin) elif select(type, rhs(_xin), `=`) <> [] then _par, _ini := selectremove(type, rhs(_xin), `=`) elif nops(rhs(_xin)) < nops(_pars)+1 then error "insufficient data for specification of initial and parameters" else _par := rhs(_xin)[-nops(_pars) .. -1]; _ini := rhs(_xin)[1 .. -nops(_pars)-1] end if; _xout := lhs(_xout); if _par <> [] then `dsolve/numeric/process_parameters`(_n, _pars, _par, _y0) end if; if _ini <> [] then `dsolve/numeric/process_initial`(_n-_ne, _ini, _y0, _pars, _vmap) end if; `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars); if _Env_smart_dsolve_numeric = true and type(_y0[0], 'numeric') and _dtbl[1][4][10] <> 1 then procname("right") := _y0[0]; procname("left") := _y0[0] end if; if _xout = "initial" then return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)] elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] else return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)], [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] end if elif _xin = "eventstop" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then return 0 end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 <= _dtbl[5-_i][4][9] then _i := 5-_i; _dtbl[4] := _i; _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) elif 100 <= _dtbl[_i][4][9] then _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) else return 0 end if elif _xin = "eventstatus" then if _nv = 0 then error "this solution has no events" end if; _i := [selectremove(proc (a) options operator, arrow; _dtbl[1][3][1][a, 7] = 1 end proc, {seq(_j, _j = 1 .. round(_dtbl[1][3][1][_nv+1, 1]))})]; return ':-enabled' = _i[1], ':-disabled' = _i[2] elif _xin = "eventclear" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then error "no events to clear" end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 < _dtbl[5-_i][4][9] then _dtbl[4] := 5-_i; _i := 5-_i end if; if _dtbl[_i][4][9] < 100 then error "no events to clear" elif _nv < _dtbl[_i][4][9]-100 then error "event error condition cannot be cleared" else _j := _dtbl[_i][4][9]-100; if irem(round(_dtbl[_i][3][1][_j, 4]), 2) = 1 then error "retriggerable events cannot be cleared" end if; _j := round(_dtbl[_i][3][1][_j, 1]); for _k to _nv do if _dtbl[_i][3][1][_k, 1] = _j then if _dtbl[_i][3][1][_k, 2] = 3 then error "range events cannot be cleared" end if; _dtbl[_i][3][1][_k, 8] := _dtbl[_i][3][1][_nv+1, 8] end if end do; _dtbl[_i][4][17] := 0; _dtbl[_i][4][9] := 0; if _dtbl[1][4][10] = 1 then if _i = 2 then try procname(procname("left")) catch: end try else try procname(procname("right")) catch: end try end if end if end if; return elif type(_xin, `=`) and member(lhs(_xin), {"eventdisable", "eventenable"}) then if _nv = 0 then error "this solution has no events" end if; if type(rhs(_xin), {('list')('posint'), ('set')('posint')}) then _i := {op(rhs(_xin))} elif type(rhs(_xin), 'posint') then _i := {rhs(_xin)} else error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; if select(proc (a) options operator, arrow; _nv < a end proc, _i) <> {} then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _k := {}; for _j to _nv do if member(round(_dtbl[1][3][1][_j, 1]), _i) then _k := `union`(_k, {_j}) end if end do; _i := _k; if lhs(_xin) = "eventdisable" then _dtbl[4] := 0; _j := [evalb(assigned(_dtbl[2]) and member(_dtbl[2][4][17], _i)), evalb(assigned(_dtbl[3]) and member(_dtbl[3][4][17], _i))]; for _k in _i do _dtbl[1][3][1][_k, 7] := 0; if assigned(_dtbl[2]) then _dtbl[2][3][1][_k, 7] := 0 end if; if assigned(_dtbl[3]) then _dtbl[3][3][1][_k, 7] := 0 end if end do; if _j[1] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[2][3][4][_k, 1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to defined init `, _dtbl[2][3][4][_k, 1]); _dtbl[2][3][1][_k, 8] := _dtbl[2][3][4][_k, 1] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to rate hysteresis init `, _dtbl[2][5][24]); _dtbl[2][3][1][_k, 8] := _dtbl[2][5][24] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to initial init `, _x0); _dtbl[2][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to fireinitial init `, _x0-1); _dtbl[2][3][1][_k, 8] := _x0-1 end if end do; _dtbl[2][4][17] := 0; _dtbl[2][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("left")) end if end if; if _j[2] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[3][3][4][_k, 2], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to defined init `, _dtbl[3][3][4][_k, 2]); _dtbl[3][3][1][_k, 8] := _dtbl[3][3][4][_k, 2] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to rate hysteresis init `, _dtbl[3][5][24]); _dtbl[3][3][1][_k, 8] := _dtbl[3][5][24] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to initial init `, _x0); _dtbl[3][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to fireinitial init `, _x0+1); _dtbl[3][3][1][_k, 8] := _x0+1 end if end do; _dtbl[3][4][17] := 0; _dtbl[3][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("right")) end if end if else for _k in _i do _dtbl[1][3][1][_k, 7] := 1 end do; _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); _dtbl[4] := 0; if _dtbl[1][4][10] = 1 then if _x0 <= procname("right") then try procname(procname("right")) catch: end try end if; if procname("left") <= _x0 then try procname(procname("left")) catch: end try end if end if end if; return elif type(_xin, `=`) and lhs(_xin) = "eventfired" then if not type(rhs(_xin), 'list') then error "'eventfired' must be specified as a list" end if; if _nv = 0 then error "this solution has no events" end if; if _dtbl[4] <> 2 and _dtbl[4] <> 3 then error "'direction' must be set prior to calling/setting 'eventfired'" end if; _i := _dtbl[4]; _val := NULL; if not assigned(_EnvEventRetriggerWarned) then _EnvEventRetriggerWarned := false end if; for _k in rhs(_xin) do if type(_k, 'integer') then _src := _k elif type(_k, 'integer' = 'anything') and type(evalf(rhs(_k)), 'numeric') then _k := lhs(_k) = evalf[max(Digits, 18)](rhs(_k)); _src := lhs(_k) else error "'eventfired' entry is not valid: %1", _k end if; if _src < 1 or round(_dtbl[1][3][1][_nv+1, 1]) < _src then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _src := {seq(`if`(_dtbl[1][3][1][_j, 1]-_src = 0., _j, NULL), _j = 1 .. _nv)}; if nops(_src) <> 1 then error "'eventfired' can only be set/queried for root-finding events and time/interval events" end if; _src := _src[1]; if _dtbl[1][3][1][_src, 2] <> 0. and _dtbl[1][3][1][_src, 2]-2. <> 0. then error "'eventfired' can only be set/queried for root-finding events and time/interval events" elif irem(round(_dtbl[1][3][1][_src, 4]), 2) = 1 then if _EnvEventRetriggerWarned = false then WARNING(`'eventfired' has no effect on events that retrigger`) end if; _EnvEventRetriggerWarned := true end if; if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then _val := _val, undefined elif type(_dtbl[_i][3][4][_src, _i-1], 'undefined') or _i = 2 and _dtbl[2][3][1][_src, 8] < _dtbl[2][3][4][_src, 1] or _i = 3 and _dtbl[3][3][4][_src, 2] < _dtbl[3][3][1][_src, 8] then _val := _val, _dtbl[_i][3][1][_src, 8] else _val := _val, _dtbl[_i][3][4][_src, _i-1] end if; if type(_k, `=`) then if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then error "cannot set event code for a rate hysteresis event" end if; userinfo(3, {'events', 'eventreset'}, `manual set event code `, _src, ` to value `, rhs(_k)); _dtbl[_i][3][1][_src, 8] := rhs(_k); _dtbl[_i][3][4][_src, _i-1] := rhs(_k) end if end do; return [_val] elif type(_xin, `=`) and lhs(_xin) = "direction" then if not member(rhs(_xin), {-1, 1, ':-left', ':-right'}) then error "'direction' must be specified as either '1' or 'right' (positive) or '-1' or 'left' (negative)" end if; _src := `if`(_dtbl[4] = 2, -1, `if`(_dtbl[4] = 3, 1, undefined)); _i := `if`(member(rhs(_xin), {1, ':-right'}), 3, 2); _dtbl[4] := _i; _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if; return _src elif _xin = "eventcount" then if _dtbl[1][3][1] = 0 or _dtbl[4] <> 2 and _dtbl[4] <> 3 then return 0 else return round(_dtbl[_dtbl[4]][3][1][_nv+1, 12]) end if else return "procname" end if end if; if _xout = _x0 then return [_x0, seq(evalf(_dtbl[1][6][_vmap[_i]]), _i = 1 .. _n-_ne)] end if; _i := `if`(_x0 <= _xout, 3, 2); if _xin = "last" and 0 < _dtbl[_i][4][9] and _dtbl[_i][4][9] < 100 then _dat := eval(_dtbl[_i], 2); _j := _dat[4][20]; return [_dat[11][_j, 0], seq(_dat[11][_j, _vmap[_i]], _i = 1 .. _n-_ne-_nd), seq(_dat[8][1][_vmap[_i]], _i = _n-_ne-_nd+1 .. _n-_ne)] end if; if not type(_dtbl[_i], 'array') then _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if end if; if _xin <> "last" then if 0 < 0 then if `dsolve/numeric/checkglobals`(op(_dtbl[1][14]), _pars, _n, _y0) then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars, _i) end if end if; if _dtbl[1][4][7] = 0 then error "parameters must be initialized before solution can be computed" end if end if; _dat := eval(_dtbl[_i], 2); _dtbl[4] := _i; try _src := `dsolve/numeric/SC/IVPrun`(_dat, _xout) catch: userinfo(2, `dsolve/debug`, print(`Exception in solnproc:`, [lastexception][2 .. -1])); error end try; if _src = 0 and 100 < _dat[4][9] then _val := _dat[3][1][_nv+1, 8] else _val := _dat[11][_dat[4][20], 0] end if; if _src <> 0 or _dat[4][9] <= 0 then _dtbl[1][5][1] := _xout else _dtbl[1][5][1] := _val end if; if _i = 3 and _val < _xout then Rounding := -infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further right of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further right of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further right of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further right of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further right of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further right of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further right of %1", evalf[8](_val) end if elif _i = 2 and _xout < _val then Rounding := infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further left of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further left of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further left of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further left of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further left of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further left of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further left of %1", evalf[8](_val) end if end if; if _EnvInFsolve = true then _dig := _dat[4][26]; _dat[4][26] := _EnvDSNumericSaveDigits; _Env_dsolve_SC_native := true; if _dat[4][25] = 1 then _i := 1; _dat[4][25] := 2 else _i := _dat[4][25] end if; _val := `dsolve/numeric/SC/IVPval`(_dat, _xout, _src); _dat[4][25] := _i; _dat[4][26] := _dig; [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] else Digits := _dat[4][26]; _val := `dsolve/numeric/SC/IVPval`(eval(_dat, 2), _xout, _src); [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] end if end proc, (2) = Array(0..0, {}), (3) = [t, x(t)], (4) = [xlow = xlow]}); _vars := _dat[3]; _pars := map(rhs, _dat[4]); _n := nops(_vars)-1; _solnproc := _dat[1]; if not type(_xout, 'numeric') then if member(x_rosenbrock, ["start", 'start', "method", 'method', "left", 'left', "right", 'right', "leftdata", "rightdata", "enginedata", "eventstop", 'eventstop', "eventclear", 'eventclear', "eventstatus", 'eventstatus', "eventcount", 'eventcount', "laxtol", 'laxtol', "numfun", 'numfun', NULL]) then _res := _solnproc(convert(x_rosenbrock, 'string')); if 1 < nops([_res]) then return _res elif type(_res, 'array') then return eval(_res, 1) elif _res <> "procname" then return _res end if elif member(x_rosenbrock, ["last", 'last', "initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x_rosenbrock, 'string'); _res := _solnproc(_xout); if _xout = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] end if elif type(_xout, `=`) and member(lhs(_xout), ["initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(lhs(x_rosenbrock), 'string') = rhs(x_rosenbrock); if type(rhs(_xout), 'list') then _res := _solnproc(_xout) else error "initial and/or parameter values must be specified in a list" end if; if lhs(_xout) = "initial" then return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] elif lhs(_xout) = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["eventdisable", 'eventdisable', "eventenable", 'eventenable', "eventfired", 'eventfired', "direction", 'direction', NULL]) then return _solnproc(convert(lhs(x_rosenbrock), 'string') = rhs(x_rosenbrock)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _vars end if; if procname <> unknown then return ('procname')(x_rosenbrock) else _ndsol; _ndsol := pointto(_dat[2][0]); return ('_ndsol')(x_rosenbrock) end if end if; try _res := _solnproc(_xout); [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] catch: error end try end proc

  

"cannot evaluate the solution further right of %1, probably a singularity", .10469398

  

 
> 

aswap := plottools[transform]((x,y) -> [y,x]):
plots[display](p, swap(p))

 
> 

 


 

Download A_Few_Ideas.mw

restart: randomize(): roll_1_dice := ra...

mmcdara 7891
October 31 2020
1 0 
restart:
randomize():
roll_1_dice := rand(1..6):
roll_3_dices := [seq(roll_1_dice(), i=1..3)]:

# a toss
toss := roll_3_dices();
# square or square root ?
map(u -> if u::odd then u^2 else sqrt(u) end if, toss);
# product
mul(toss);


But perhaps it would be more beneficial to you if you tried to make your own?
Don't you think?

Your problem is not completely defined. ...

mmcdara 7891
October 31 2020
1 1

Your problem is not completely defined.
There exist two ways to understand it:
 

> 

restart:
> 

# First question : are your dices distinct or not ?
> 

randomize():
> 

roll_1_dice := rand(1..6):
> 

roll_2_dices := [roll_1_dice(), roll_1_dice()]

[4, 3]

 (1)
> 

# First case : distinct dice

# Sequence of all outcomes
# There are many ways to this, for instance


with(combinat, cartprod):
T := cartprod([[$1..6], [$1..6]]):
all_outcomes := NULL:
while not T[finished] do
  all_outcomes := all_outcomes, T[nextvalue]()
end do:
all_outcomes_D := [all_outcomes]

[[1, 1], [1, 2], [1, 3], [1, 4], [1, 5], [1, 6], [2, 1], [2, 2], [2, 3], [2, 4], [2, 5], [2, 6], [3, 1], [3, 2], [3, 3], [3, 4], [3, 5], [3, 6], [4, 1], [4, 2], [4, 3], [4, 4], [4, 5], [4, 6], [5, 1], [5, 2], [5, 3], [5, 4], [5, 5], [5, 6], [6, 1], [6, 2], [6, 3], [6, 4], [6, 5], [6, 6]]

 (2)
> 

# another way

M := Matrix(6, 6):
all_outcomes_D := [indices(M)]

[[1, 1], [2, 1], [3, 1], [4, 1], [5, 1], [6, 1], [1, 2], [2, 2], [3, 2], [4, 2], [5, 2], [6, 2], [1, 3], [2, 3], [3, 3], [4, 3], [5, 3], [6, 3], [1, 4], [2, 4], [3, 4], [4, 4], [5, 4], [6, 4], [1, 5], [2, 5], [3, 5], [4, 5], [5, 5], [6, 5], [1, 6], [2, 6], [3, 6], [4, 6], [5, 6], [6, 6]]

 (3)
> 

# still another one

all_outcomes_D := [seq(seq([i, j], j=1..6), i=1..6)]

[[1, 1], [1, 2], [1, 3], [1, 4], [1, 5], [1, 6], [2, 1], [2, 2], [2, 3], [2, 4], [2, 5], [2, 6], [3, 1], [3, 2], [3, 3], [3, 4], [3, 5], [3, 6], [4, 1], [4, 2], [4, 3], [4, 4], [4, 5], [4, 6], [5, 1], [5, 2], [5, 3], [5, 4], [5, 5], [5, 6], [6, 1], [6, 2], [6, 3], [6, 4], [6, 5], [6, 6]]

 (4)
> 

# Randomly pick an outcome
# Here again a solution among many others

combinat:-randperm(all_outcomes_D)[1]

[2, 1]

 (5)
> 

#  another solution

pick_D := rand(1..numelems(all_outcomes_D)): # equivalently rand(1..36):
all_outcomes_D[pick_D()]

[5, 3]

 (6)
> 

# Second case : undistinct dice

all_outcomes_U := [seq(seq([i, j], j=i..6), i=1..6)]

[[1, 1], [1, 2], [1, 3], [1, 4], [1, 5], [1, 6], [2, 2], [2, 3], [2, 4], [2, 5], [2, 6], [3, 3], [3, 4], [3, 5], [3, 6], [4, 4], [4, 5], [4, 6], [5, 5], [5, 6], [6, 6]]

 (7)
> 

#  Randomly pick an outcome

pick_U := rand(1..numelems(all_outcomes_U)):
all_outcomes_U[pick_U()]

[3, 4]

 (8)
> 

# test (distinct dice)

toss   := mul(roll_2_dices());
picked := mul(all_outcomes_D[pick_D()]);

if   toss < picked then "less"
elif toss = picked then "equal"
else                    "greater"
end if;

12

  

16

  

"less"

 (9)
> 

 


 

Download A_solution.mw

@9009134 For information here is th...

mmcdara 7891
September 27 2020
0 14

@9009134

 

For information here is the result the statistical software R produces (fit of CP/R by a 4th degree polynomial in T)

t <- c(298., 300., 400., 500., 600., 700., 800., 900., 1000.)
cp <- c(54.187, 54.518 , 70.421 , 84.683 , 96.982 , 107.323, 116.042, 123.801 , 131.594)
cpoverR <- cp / 1.987
lm(cpoverR ~ 1+t+I(t^2)+I(t^3)+I(t^4))  # linear fit

(Intercept)            t                    I(t^2)           I(t^3)            I(t^4)  
 2.579e+00    7.263e-02    7.326e-05   -1.512e-07    6.895e-11  

Which means the model R finds is 2.579 + 7.263e-02*t + 7.326e-05*t^2 -1.512e-07*t^3 + 6.895e-11*t^4
Practically identical to the model I found "by hand".
 


A way to force LinearFit is to operate on SCALED data.
Here is an illustration for the model CP/R = f(T)
Now LinearFit and the "by hand" method (such as R) give the same result.

A remark for people from the development team that could read your thread.
I consider that a clever implementation of LinearFit should account for ill condionned problems.
Solutions are well known 

  • normalize the data (as I did here)
  • use methods less sensitive to bad condionning (Givens-Householder for instance)
  • use othogonal polynomials and express the result in canonical basis
     
> 

restart:
> 

with(Statistics):
> 

T := [298., 300., 400., 500., 600., 700., 800., 900., 1000.];
> 

CP := [54.187, 54.518 , 70.421 , 84.683 , 96.982 , 107.323, 116.042, 123.801 , 131.594];

[298., 300., 400., 500., 600., 700., 800., 900., 1000.]

  

[54.187, 54.518, 70.421, 84.683, 96.982, 107.323, 116.042, 123.801, 131.594]

 (1)
> 

R := 1.987;

1.987

 (2)
> 

# Let's try to scale the data (case of T and CP only)

Tscaled  := Scale(T):
CPscaled := Scale(CP/~R):

Tscaled, CPscaled;

`Cp/R`    := theta -> theta^4*a5+theta^3*a4+theta^2*a3+theta*a2+a1;

M1scaled := LinearFit(`Cp/R`(theta), Tscaled, CPscaled, theta);

Matrix(9, 1, {(1, 1) = HFloat(-1.2154612619718141), (2, 1) = HFloat(-1.207691978337051), (3, 1) = HFloat(-0.8192277965989003), (4, 1) = HFloat(-0.43076361486074943), (5, 1) = HFloat(-0.04229943312259854), (6, 1) = HFloat(0.34616474861555213), (7, 1) = HFloat(0.7346289303537028), (8, 1) = HFloat(1.1230931120918535), (9, 1) = HFloat(1.5115572938300046)}), Matrix(9, 1, {(1, 1) = HFloat(-1.3453854645769592), (2, 1) = HFloat(-1.333995104145723), (3, 1) = HFloat(-0.7867416300584438), (4, 1) = HFloat(-0.2959581904637911), (5, 1) = HFloat(0.12727456316287933), (6, 1) = HFloat(0.48312868952923926), (7, 1) = HFloat(0.7831666091417269), (8, 1) = HFloat(1.0501690435465787), (9, 1) = HFloat(1.318341483864494)})

  

proc (theta) options operator, arrow; theta^4*a5+theta^3*a4+theta^2*a3+theta*a2+a1 end proc

  

HFloat(0.16917822809413155)+HFloat(0.9812133698557858)*theta-HFloat(0.22397858937737267)*theta^2+HFloat(0.020173883967869902)*theta^3+HFloat(0.020702898675047213)*theta^4

 (3)
> 

# Back to the original variables

ScaledCPoverR := eval(M1scaled, [theta = (theta - Mean(T)) / StandardDeviation(T)]);

CPoverR := ScaledCPoverR * StandardDeviation(CP/~R) + Mean(CP/~R);

-HFloat(2.159324034212451)+HFloat(0.003811662488315613)*theta-HFloat(0.22397858937737267)*(HFloat(0.003884641817381508)*theta-HFloat(2.3730845235515035))^2+HFloat(0.020173883967869902)*(HFloat(0.003884641817381508)*theta-HFloat(2.3730845235515035))^3+HFloat(0.020702898675047213)*(HFloat(0.003884641817381508)*theta-HFloat(2.3730845235515035))^4

  

HFloat(15.366996346473783)+HFloat(0.05574515184793432)*theta-HFloat(3.275662657384246)*(HFloat(0.003884641817381508)*theta-HFloat(2.3730845235515035))^2+HFloat(0.2950408722175395)*(HFloat(0.003884641817381508)*theta-HFloat(2.3730845235515035))^3+HFloat(0.3027776551231061)*(HFloat(0.003884641817381508)*theta-HFloat(2.3730845235515035))^4

 (4)
> 

M1unscaled:= expand(%);

HFloat(2.5793789487674212)+HFloat(0.07262797433002269)*theta+HFloat(7.325618169228684e-5)*theta^2-HFloat(1.511850442093991e-7)*theta^3+HFloat(6.89489640013437e-11)*theta^4

 (5)
> 

span            := unapply(subs([seq(a||i=1, i=1..5)], [op(`Cp/R`(theta))]), theta):
DesignMatrix    := convert(span~(T), Matrix):
D2Matrix        := (DesignMatrix^+ . DesignMatrix)^(-1):
ByHandSolution  := Vector(5, i -> a||(6-i))
                   =
                   D2Matrix . DesignMatrix^+ . convert(CP/~R, Vector):

M1ByHand := convert(span(theta), Vector[row]) . rhs(ByHandSolution);

HFloat(2.579378948829685)+HFloat(6.894896400274164e-11)*theta^4-HFloat(1.5118504421633464e-7)*theta^3+HFloat(7.325618168985484e-5)*theta^2+HFloat(0.07262797432690293)*theta

 (6)
> 

plots:-display(
  ScatterPlot(T, CP/~R),
  plot(M1unscaled, theta=min(T)..max(T), color=blue, legend='LinearFit'),
  plot(M1ByHand, theta=min(T)..max(T), color=red, legend='ByHand'),
  gridlines=true
);

 
> 

plot(M1unscaled-M1ByHand, theta=min(T)..max(T), title="difference between models", gridlines=true)

 
> 

 


 

Download FITTING_after_scaling.mw

A solution using the ordinary least squa...

mmcdara 7891
September 27 2020
0 0

Here are fittings obtained the least squares method.

You will see that the Maple procedure LinearFit suffers a lot of problems :

  1. It retuns a warning saying the (information) matrix is not of full rank.
    I show it is the case for Digits=10 (default value), but that this same matrix is of full rank for a large enough value of Digits (here 40). In my opinion this denotes a weakness of LinearFit
  2. For the model S/T=f(T) (written this way, but a clever way would be to fit S=g(T) and next write S/T = g(T)/T...) LinearFit returns a completely absurd error (Error, (in Statistics:-LinearFit) the model, `S/T)`(theta), is not linear in the parameters, [])

So I completed the use of LInearFit with a "By hand" construction of the models (an extremely easy thing to do if you have basic knowledge about the ordinary least squares method).
The "rank problem" can be easily avoided and no "error" (as expected) occur when proceeding this way.

In my opinion I consider LinearFit suffers some weaknesses. 

I didn't try to use NonlinearFit instead of linearFit as a workaround because it is just stupid given all your models are linear.

Download FITTING.mw

Now, a little advice about the "good practices in linear regression".
A good way to fit statistical models is to begin by looking at the data and infer some simple model (for instance a low degree polynomial model).
If you do thys for the CP/R=f(T) model you will see it's mainly linear with only a slight discrepancy that indicates adding a quadratic term could be sufficient.
So, try to fit CP/R to the model a0+a1*T and if the result doesn't suit you fit the model a0+a1*T*a2*T^2 (and augment it if you want.

Of course (assuming no warnong or error occurs) the higher the number of parametrers (= 1+degree of the model), the lower the residual sum of squares (RSS), providing this number doesn't exceed the number of observations (9 in your case).
Let's suppose you fit polynomial models of degree d and denote by RSS(d) the value of RSS for the polynomial model of degree d ; then RSS(0) > RSS(1) > ... > RSS(8) = 0 (with your 9 data).
But be extremely careful, RSS is NOT the good criteria to use to select the best model among these d-1 models.

...

mmcdara 7891
September 19 2020
1 0

No other solution than the trivial one if sin(y) <> 0

PS : Of course, if y=k*Pi  (k in Z) , then any function which respects the boundary conditions is solution. for instance GenFunc(x) = (1-cos(x*2*Pi/a))
 

``

> 

restart:
> 

w := sin(y)*GenFunc(x)

sin(y)*GenFunc(x)

 (1)
> 

# given all the dij are identical, Nx=Ny and Nxy=0 we have

EDO := D__11*diff(w, x$4)+6*D__11*diff(w, [x$2, y$2])+D__11*(diff(w, y$4))+N__x*(diff(w, x$2))+N__x*(diff(w, y$2))

D__11*sin(y)*(diff(diff(diff(diff(GenFunc(x), x), x), x), x))-6*D__11*sin(y)*(diff(diff(GenFunc(x), x), x))+D__11*sin(y)*GenFunc(x)+N__x*sin(y)*(diff(diff(GenFunc(x), x), x))-N__x*sin(y)*GenFunc(x)

 (2)
> 

EDO := simplify(EDO)

sin(y)*(D__11*(diff(diff(diff(diff(GenFunc(x), x), x), x), x))-6*D__11*(diff(diff(GenFunc(x), x), x))+D__11*GenFunc(x)+N__x*(diff(diff(GenFunc(x), x), x))-N__x*GenFunc(x))

 (3)
> 

EDO := op(2, EDO);

D__11*(diff(diff(diff(diff(GenFunc(x), x), x), x), x))-6*D__11*(diff(diff(GenFunc(x), x), x))+D__11*GenFunc(x)+N__x*(diff(diff(GenFunc(x), x), x))-N__x*GenFunc(x)

 (4)
> 

# To see what happens, set boundary conditions at x=0 alone.

SOL := unapply(rhs(dsolve({EDO=0, GenFunc(0) = 0, D(GenFunc)(0) = 0})), x);
 

proc (x) options operator, arrow; (-(1/2)*((D__11*(6*D__11-N__x-(32*D__11^2-8*D__11*N__x+N__x^2)^(1/2)))^(1/2)+(D__11*(6*D__11-N__x+(32*D__11^2-8*D__11*N__x+N__x^2)^(1/2)))^(1/2))*_C3/(D__11*(6*D__11-N__x-(32*D__11^2-8*D__11*N__x+N__x^2)^(1/2)))^(1/2)-(1/2)*((D__11*(6*D__11-N__x-(32*D__11^2-8*D__11*N__x+N__x^2)^(1/2)))^(1/2)-(D__11*(6*D__11-N__x+(32*D__11^2-8*D__11*N__x+N__x^2)^(1/2)))^(1/2))*_C4/(D__11*(6*D__11-N__x-(32*D__11^2-8*D__11*N__x+N__x^2)^(1/2)))^(1/2))*exp(-(1/2)*2^(1/2)*(D__11*(6*D__11-N__x-(32*D__11^2-8*D__11*N__x+N__x^2)^(1/2)))^(1/2)*x/D__11)+(-(1/2)*((D__11*(6*D__11-N__x-(32*D__11^2-8*D__11*N__x+N__x^2)^(1/2)))^(1/2)-(D__11*(6*D__11-N__x+(32*D__11^2-8*D__11*N__x+N__x^2)^(1/2)))^(1/2))*_C3/(D__11*(6*D__11-N__x-(32*D__11^2-8*D__11*N__x+N__x^2)^(1/2)))^(1/2)-(1/2)*((D__11*(6*D__11-N__x-(32*D__11^2-8*D__11*N__x+N__x^2)^(1/2)))^(1/2)+(D__11*(6*D__11-N__x+(32*D__11^2-8*D__11*N__x+N__x^2)^(1/2)))^(1/2))*_C4/(D__11*(6*D__11-N__x-(32*D__11^2-8*D__11*N__x+N__x^2)^(1/2)))^(1/2))*exp((1/2)*2^(1/2)*(D__11*(6*D__11-N__x-(32*D__11^2-8*D__11*N__x+N__x^2)^(1/2)))^(1/2)*x/D__11)+_C3*exp(-(1/2)*2^(1/2)*(D__11*(6*D__11-N__x+(32*D__11^2-8*D__11*N__x+N__x^2)^(1/2)))^(1/2)*x/D__11)+_C4*exp((1/2)*2^(1/2)*(D__11*(6*D__11-N__x+(32*D__11^2-8*D__11*N__x+N__x^2)^(1/2)))^(1/2)*x/D__11) end proc

 (5)
> 

# Obviously Genfunc(x) is not identically null unless _C3=_C4=0

evalf(eval(SOL(0.01), {D__11=10000, N__x=1000}));
evalf(eval(SOL(0.02), {D__11=10000, N__x=1000}));

0.277112e-3*_C3+0.281574e-3*_C4

  

0.1099684e-2*_C3+0.1135393e-2*_C4

 (6)
> 

# let's find _C3 and _C4 such that SOL(a)=D(SOL)(a)=0 :

DSOL := unapply(diff(SOL(x), x), x):
solve({SOL(a)=0, DSOL(a)=0}, {_C3, _C4})

{_C3 = 0, _C4 = 0}

 (7)
> 

# thus the only solution is GenFunc(0)=0 ... QED

0

 (8)
> 

 


 

Download N0_Other_Solution.mw

Here is something wich looks like what y...

mmcdara 7891
September 19 2020
2 2

Here is something wich looks like what you want.
Colors are randomly choosen, so the result is not always pretty (there is a randomize() command at the top of the worksheet to generate a random seed which will enable to obtain different colourings). 
You can of course define your own colors.

To obtain something closer to the picture tou present it would be necessary to rescale the radii of the vertices of each graph, a thing I didn't do here.

Left : graph1 ; middle : graph 2 ; right carteian product of g1 and g2.
Cartesian_Product_of_Graphs.mw

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