mmcdara

7891 Reputation

22 Badges

9 years, 52 days

MaplePrimes Activity


These are replies submitted by mmcdara

@Kitonum 

I'm so disappointed not to have thought of it myself, it was such a simple idea.
Shame on me and bravo to you.

@ashy 

I use above the terme sensitivity without any definition.
For a reference on the subject I advice you to refer to the link I provide at the top of the attached file.
Note that Wikipedia already provides simple definitions.

As you assess the sensitivity of Ro to some parameter p by using partial derivatives, this means you are doing what is called local sensitivity (by opposition to global sensitivity [see these terms on Wikipedia for instance]).

In local sensitivity the variations of the parameters are aimed to be infinetisimal (simply because your computations are based on partial derivatives). The signs of the partial derivatives give you an indication in the direction of variation of Ro at the reference point you use, in the "p direction" (does Ro increases as p increases too, or does Ro decreases as p increases?) and the absolute values of these derivatives give you an indication of the magnitude of variation of Ro with respect to p.
I think we both agree on this point.

The local variation of some function f(p, q) at point [p=a, q=b] can be obtained by usng a simple Taylor exansion:

df := unapply( mtaylor(f(p, q), [p=a, q=b], 2), (a, b))
    df := (a, b) -> f(a, b) + D[1](f)(a, b) (p - a) + D[2](f)(a, b) (q - b)

The part of df(a, b) which comes from, let's say p, is equal to D[1](f)(a, b) (p - a) where p - a is assumed to be small:

  • The first point is that you do not write something like this but you consider D[1](f)(a, b)*a  instead.
    What is the motivation in multiplying by a?
  • The second point is that you introduce a second "correction" to write finally D[1](f)(a, b)*a / f(a, b).
    Why do you divide by f(a, b)  (you could have divide by any other non null constant the same way)?
    What is perturbing is that doing so you do not compute the partial derivative of f but those of log(f):why this unnecessary complication?

When I say "Maybe they are of common use in epidemiology?" I mean

  • Is it usual in epidemiology to assess the sensitivities of some function f by computing, not its partial derivatives, but those of log(f) (which would justify you writing diff(Ro, p)/Ro = diff(log(Ro), p)*p?

Here is a recent article about Sansitivity analysis of SIR-type models, read chapter B. Sensitivity Analysis Methods, more specifically paragraph 1) Local Sensitivity Analysis : "A local sensitivity index is mathematically the partial derivative of model outputs concerning an individual parameter" (note this definition does not conform the definition you will find following the link given in the attached file as the the term index should have been replaced by coefficient), so, in this specific paper that I'm not going to pretend is an undisputable reference, D[1](f)(a, b) is a "sensitivity index" but is D[1](f)(a, b)*a / f(a, b) not.

So, indeed, the definition you use is debatable.

discussion.mw

@salim-barzani 

As I told you here it's better fo work with the initial simple expression as long as you can and, once you have got to the results you are looking for to make the substitutions you want, rather tha complexifying the initial expression by making these substitutions first and then trying to get the results.

You must understand that M is an amazingly complex function of x (or y, or t) depending on 16 parameters and that you have absolutly no chance to find combinations of these parameters in order that this function has no singularities.

One can show that singularities are the values of x that verify 

Sum(Ri*exp(Ai*x+Bi), i=1..15) = 0

where Ri, Ai, Bi both depend on the 16 parameters (and on y and t).
What you want is finding those 16 parameters in order that there exist no x value which verify the above equation.
This is complete illusion.

So either you provide a more synthetic expression of M (the native one), or, if you are lucky enough, you must adopt a brute force approach where you randomly chose parameter values and look if they lead to nonsingular solutions.

The attached file explores that last option

MAPLEPRIMES ISSUE : < select file  not responding>
I can't attach the file, I will try again ia a few minutes.


Here it is, file attached (Does Mapleprimes fully recover after its breakout?)
Meanwhile here is an exercpt

# With integer rand values in between -1 and +1

found := false:

while not found do
  try
    data := params =~ [seq(rand(-1 .. 1)(), k=1..numelems(params))]:
    DATA := [data[], y=0, t=0];

    ux := eval(M, DATA);

    fd := fdiscont(ux, x=-10..10):
    if fd = [] then found := true end if:
  catch:
  end try;
end do:

data;
   [a = 0, alpha = 1, b = 1, beta = -1, a[1] = -1, a[2] = 0, 

     a[3] = -1, b[1] = 1, b[2] = -1, b[3] = 1, eta[3] = 0, 

     eta[4] = 0, k[3] = 1, k[4] = 0, l[3] = 0, l[4] = -1]

ux;
(2 (-2 I exp(-I x + (-1 + I)) + 2 I exp(I x + (-1 - I))

   + (-1 - I) exp(I x + (-1 - I) + x) + 2 exp(-2 + x)

   + (-1 + I) exp(-I x + (-1 + I) + x) + 2 exp(x)))/(2

   + 2 exp(-I x + (-1 + I)) + 8 exp(-2) + 2 exp(I x + (-1 - I))

   - exp(I x + (-1 - I) + x) + 2 exp(-2 + x)

   - exp(-I x + (-1 + I) + x) + 2 exp(x))


DocumentTools:-Tabulate( 
  [    
    plot(ux),  
    plot3d(
      eval(M, [data[], t=0])
      , y = -10 .. 10
      , x = -10 .. 10
      , shading = ZHUE
      , grid = [500, 500]
      , lightmodel = light4
      , style = surfacecontour
    ) 
  ]   
  , width=50 
):


singular-interval_mmcdara(1).mw
 

Examples:
A solution with integer valued parameters in the set {-1, 0, 1}





A solution with float valued parameters in the interval [-1, 1]

@salim-barzani 

M depends on those 16 parameters (x, y, t put aside)

{a, alpha, b, beta, t, x, y, a[1], a[2], a[3], b[1], b[2], b[3], eta[3], eta[4], k[3], k[4], l[3], l[4]}

and you want to find values for them in order that M has no singularities, am I right?

@Alfred_F 

I never said the serie was not convergent: I only said it was not summable.

The classical example of a convergent serie which is not summable is sum((-1)n / n, n=1..+oo): as its limit is -ln(2), splitting this serie as sum((-1)2*n-1 / (2*n-1), n=1..+oo) + sum((-1)2*n / (2*n), n=1..+oo) converges to -ln(2)/2

So we must be fery careful in the limit we (here Maple) get.

Before going any further, I have to say that I'm not a specialist, and that what I'm about to say may be incorrect.

Let A the set of couples (m, n) where both m and n are strictly positive integers such that m <> n.
Then you can separate A as A> + A<  where (m, n) belongs to A> if m > n and  (m, n) belongs to A< if m < n.
Rewritting S = sum(1/(m^2-n^2), (m, n) in A) as S = sum(1/(m^2-n^2), (m, n) in A>) + sum(1/(m^2-n^2), (m, n) in A<) gives (immediate result) S = 0.
So the limit os the serie depends on the way the summation is done
Thus depending on the way the summation is done the limits are diferent.

restart:

s0 := Sum(Sum(1/(m^2-n^2), m=1..+infinity), n=1..+infinity)

Sum(Sum(1/(m^2-n^2), m = 1 .. infinity), n = 1 .. infinity)

(1)

s1 := Sum(Sum(1/(m^2-n^2), m=1..n-1), n=1..+infinity) + Sum(Sum(1/(m^2-n^2), m=n+1..+infinity), n=1..+infinity)

Sum(Sum(1/(m^2-n^2), m = 1 .. n-1), n = 1 .. infinity)+Sum(Sum(1/(m^2-n^2), m = n+1 .. infinity), n = 1 .. infinity)

(2)

s2 := value(s1)

sum(-(1/2)*Psi(2*n)/n-(1/2)*gamma/n+(1/2)*Psi(n+1)/n-(1/2)*Psi(n)/n, n = 1 .. infinity)+sum((1/2)*Psi(2*n+1)/n+(1/2)*gamma/n, n = 1 .. infinity)

(3)

combine(s2)

(1/8)*Pi^2

(4)

t11 := Sum(Sum(1/(m^2-n^2), m=1..n-1), n=1..+infinity);  # sweeps L by row and, for each rox, by column
t12 := Sum(Sum(1/(m^2-n^2), n=1..m-1), m=1..+infinity);  # sweeps U by column and, for each column, by row

Sum(Sum(1/(m^2-n^2), m = 1 .. n-1), n = 1 .. infinity)

 

Sum(Sum(1/(m^2-n^2), n = 1 .. m-1), m = 1 .. infinity)

(5)

# Another way to do the summation

t21 := value(t11);
t22 := value(t12);

combine(subs(n=k, %%) + subs(m=k, %));

sum(-(1/2)*Psi(2*n)/n-(1/2)*gamma/n+(1/2)*Psi(n+1)/n-(1/2)*Psi(n)/n, n = 1 .. infinity)

 

sum((1/2)*Psi(2*m)/m+(1/2)*gamma/m-(1/2)*Psi(m+1)/m+(1/2)*Psi(m)/m, m = 1 .. infinity)

 

0

(6)

# Imagine an infinite matrix M with null diagonal.
#
# The first summation strategy (s1) sweeps M row by row and, for each row, column by column.
#
# The summation summation strategy:
#       (t11) sweeps M row by row and, for each row, column by column from column 1 to the diagonal of M
#       (t12) sweeps M column by column and, for each column, row by row from row 1 to the diagonal of M
#
# As 1/(m^2-n^2) = -1/(n^2-m^2), it is obvious this second strategy gives 0 as limit

# WHY THE FIRST STRATEGY (which I used in my answer as dharr did) WOULD BE THE ONE TO USE?

 ``

Download summability.mw

The limit pi2/8 is nothing but some possible limit value of S for a given arbitrary choice of a summation strategy.

@dharr @acer @Alfred_F 

The serie is not summable... which raises some questions concerning  Maple's Pi^2/8 result.

Here is a demonstration (in French): see Exercice 21, p 3 (Exercise 21) and its correction (follow the link).
Translations:

  • Justifier : Prove that
  • En déduire : Deduce that
  • Qu'en déduire? : What conclusion can we draw?

Translations (from correction)

  • La série converge compte tenu des critères usuels : Convergence comes from clasical criteria
  • Par telescopage : Using a telescoping sum
  • De plus : More of this
  • Donc : Thus
  • Puis : Next
  • Cependant : But
  • On en déduit... : We deduce that the family...is not summable.

Observe that the first question in Exercice 21 corresponds to @acer's inner sum comment

@Alfred_F 

What will happen when to the condition m <> n when m = +oo and n=+oo but m=n-p for any strictly posiive defined integer p?

restart

t := (m, n) -> piecewise(m = n, 0, 1/(m^2-n^2))

proc (m, n) options operator, arrow; piecewise(m = n, 0, 1/(m^2-n^2)) end proc

(1)

S := M -> Sum(Sum(t(m, n), m=1..M), n=1..M);

value(S(10));   # can be proved by simple symmetry argument
value(S(100));  # same thing

proc (M) options operator, arrow; Sum(Sum(t(m, n), m = 1 .. M), n = 1 .. M) end proc

 

0

 

0

(2)

# But, as N -> infinity:

N := +infinity:
t(N, N);

N-1;
is(N=N-1);
t(N, N-1)

0

 

infinity

 

true

 

0

(3)

# The problem with infinite sums and a condition like m <> n (or m=n) is that they contain
# terms like t(+oo, +oo - some_definite_integer) and that this means that wou want to
# compare infinities.
#
# So there it has no sense to use a condition such as m=n when m and n are infinite


Download Infinities.mw

@nm 

I use to use notepad++ at the office and you're right, it is a very good tool especially when you use syntax coloring and code versionning.
Unfortunately when I'm at home I use an iMac and I wasn't capable to find any "true" alternative to notepad++ for Mac OSX: are you aware os some notepad++-like for Mac OSX?

TIA

@Arya-S-AA 

I'm happy for you


want to find 7 parameters which verify 252 equations ???

numelems(eqs); 
                              252
numelems(vars);
                               7

Maybe you shoulg think twice about that, unless your et of equations has a very specific structure (for instance 245 equations are linear combinations of the remaining 7 equations, for instance)

Did you give a look to the IterativeMaps package?

What does your phrase "So we want to find a substitution that removes the time dependence from u. One way is to find the maximum and see how it moves. Here, the first solution gives what we want." mean precisely?

What kind of transformation are you talking about?

It is obvious that eq17 also writes

u(x, y, z, t) = 2*diff(L(x, y, z, t)), x$2);

# with 

L(x, y, z, t) = log(f, x, y, z, t)

So if you don't want u to depend on t you have to take f  to be for instance

f(x, y, z, t) = A(x, y, z)*B(t)

Is that what you need?

In such case have a look to  hfz_mmcdara.mw

If not forget it.


For_your_information.mw shows it worked well with Maple 2015: likely a code regression

@Kitonum 

"The variant ... more universal"

I agree but it wasn't the question.

BTW an answer coser to what the OP asked would have been

A := n -> Matrix(n, (i, j) -> a_||i||j);
 #or
A := n -> Matrix(n, (i, j) -> cat(a_, i, j))

Have a good day

@Kitonum 

Why didn't you write simply

A:=n->Matrix(n, symbol=a):

instead of 

A:=n->Matrix(n,(i,j)->a[i,j]):

?

Indeed the two definitions below seem equivalent

A:=n->Matrix(n,(i,j)->a[i,j]):
lprint(A(2))
Matrix(2, 2, {(1, 1) = a[1, 1], (1, 2) = a[1, 2], (2, 1) = a[2, 1], (2, 2) = a[2, 2]}, datatype = anything, storage = rectangular, order = Fortran_order, shape = [])

B:=n->Matrix(n, symbol=a):
lprint(B(2))
Matrix(2, 2, {(1, 1) = a[1, 1], (1, 2) = a[1, 2], (2, 1) = a[2, 1], (2, 2) = a[2, 2]}, datatype = anything, storage = rectangular, order = Fortran_order, shape = [])
4 5 6 7 8 9 10 Last Page 6 of 154