restart:

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

# h = max(f, g)

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

# example 1 (Katatonia's)

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

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

# a lengthy output not represented

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

# integration for Katatonia's example

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

# example continued

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

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

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

# example 2

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

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

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

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

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

restart

with(GraphTheory):

with(SpecialGraphs):

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

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

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

restart:

with(Statistics):

#choose a vector of 4 strictly positive values

alpha := [1$4]:

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

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

# compute N row vectors whose sum is 1

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

# verify

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

restart:

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

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

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

# 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;

# 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);
 

# Next eval d2a at LS and RS

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

# 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=%);

# note 1: this command is useful to capture discontinuities

discont(beta, h)

# 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);

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)

EgyptianFraction(5/121)

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)

BruteForce(2)

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

restart:

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

plots[interactive]();

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

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

# 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]();

# 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]);

# 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

# 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) ]

# 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)

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

s_values := F~(t_values)

# 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)^+

# 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

# 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)

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

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);

# 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);

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

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

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

# 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)));

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

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)])

restart:

interface(version)

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:

# 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):

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

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

# 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):
 

# 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):

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

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

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)

 

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)

# let's check this

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

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

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

# redefine f(n) correctly:

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

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

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

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

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

# 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) ]

# 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))));

# 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

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)

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