4500 Reputation

17 Badges

6 years, 289 days

MaplePrimes Activity

These are replies submitted by mmcdara


Here is an edited version of the worksheet I already sent you.

In particular you will see how to use Optimization:-Minimize for future needs.
The main problem comes from the fact that K is to be an integer and that NLPSolve doesn't handle integer variables.
So a strategy could be to find the minimum over (A, M) for any integervalue of K in some range 0...Kmax (see the file above).

Watchout, Optimization:-Minimize is extremely slow and I think that using it requires some preliminary simplifications of the function to minimize (if possible).

The error message seems quite clear (Maple 2015.2)

Error, (in pdsolve/numeric/process_PDEs) can only numerically solve PDE with two independent variables, got {r, t, z}

Your pde contains 3 independent variables t, r, an z and Maple can't handle this.

REMARK: event the stationnary problem (a Laplace equation in axi-symetric cylindrical coordinates) cannot be handled by Maple (2015.2) 
Error, (in pdsolve/numeric) unable to handle elliptic PDEs


And I've deleted mine.

Good point for seeing that it was volts.
I vote up


Maybe the notation A:=B+1; is a little but confusing (I didn't even know that it would give the expected result) ?
So write it this way

local I:
B := <0,-2,1;-1,0,0;-2,0,0>:
A := I+B:

If A = I +B then
A^n = I + n*B + n*(n-1)/2*B^2 + n(n-1)*(n-2)/6 + .... B^n
(apply to a sum of matrices what you do for a "simple" sum)

Observe that binomial(n, 1) = 1, binomial(n, 2) = n*(n-1)/2 and so on.

Now @vv showed that B^3 = 0 ==> B^4 = B.B^3 = 0 ... ==> B^p = 0 for each p > 2.

Thus   A^n = I + n*B + n*(n-1)/2*B^2 + n(n-1)*(n-2)/6 + .... B^n  =  I + n*B + n*(n-1)/2*B^2

'A'^n = I + n*B + binomial(n,2)*B^2


@Thomas Richard 

I didn't know about this function.
I'll have a look at it as soon as I get back to my office.
Anyway, thanks for the tip.

@Thomas Richard 

Here is an example when "capturing" this information could be useful.
Suppose you have a file which contains the definition of a module named foo  and a command 

savelib(foo, "foo.mla"); 

You execute the worksheet at date T.
How can you be 100% sure that the foo.mla you use at date T' > T has been created from the that existed at time T?
To  answer this question I decided years ago to "programatically capture" the name of the worksheet within this procedure, its modifaction (through FileTools) and the name of the operator (ssystem("whoami")) and augment the archive foo.mla  with the informations.

 Taches  := ssystem("tasklist /V /FO ""LIST"" "):
 MaTache := StringTools:-StringSplit(Taches[2], "javaw.exe")[2]:
 MaTache := StringTools:-StringSplit(%, "Maple 2019")[1]:
 MaTache := StringTools:-Squeeze(MaTache):
 MaTache := StringTools:-StringSplit(MaTache, "fenˆtre: ")[2]:
 MyWS    := StringTools:-StringSplit(MaTache, "*")[1]:
 MyWS    := StringTools:-SubstituteAll(MyWS, "\\", "/"):
 WSepoch := FileTools[ModificationTime](MyWS):

 Seed  := randomize():
 GenID := StringTools:-Random(20, "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"):
 INFO := table([
     "CreationDate"       = StringTools[FormatTime]("%Y.%m.%d_%H.%M.%S"),
     "WorksheetName"      = MyWS,
     "WorksheetEpochTime" = WSepoch,
     "Author"             = op(-1, ssystem("whoami"))

 LibraryTools:-Save(INFO, "foo.mla"):

By the time I couldn't find any versioning tools that were secure enough to guarantee the link between and foo.mla.

The weak link in my strategy is that you heve to save before executing it and thus creating foo.mla.


Very nice solution, I vote up!

I was working on a solution based upon RegSub/RegSubs.
Both of these functions seem very powerful but I find their help pages somewhat obscure.
Could you give me some hints on how they work?



You're right, the qolution is OS-dependent.
For instance, for years I use the trick I explain in my answer to @mjc .

A priori, an equivalent strategy could be based on this command if you use a Linux machine, but I never checked this.

xprop -id $(xprop -root -f _NET_ACTIVE_WINDOW 0x " \$0\\n" _NET_ACTIVE_WINDOW | awk "{print \$2}")

The thing is it doesn't work on Mac OSX (_NET_ACTIVE_WINDOW is unedfined) and I never found out how to reproduce with Mac OSX what I do with Windows (no clear equivalent to tasklist).

Maybe this can give you some ideas if you want to go further


Thanks for your reply.
At least this shows that there is no version issue.


"and there is probably not much you can do about it!"
You're wrong, simply the resolution of hyperbolic pde requires a minimum of precautions.
I invite you to consult my answer.

Nevertheless, I do acknowledge that there may be a version problem (2015 = 2020 <> yours).

@Athar Kharal 

See my answer from my home login @mmcdara.
When I was at the office (@sand15) I used Maple 2020 and it worked just fine (like Maple 2015 if I'm not mistaken).


It's a matter of point of view and probably of usage. and historical conventions.

Take the Boyle-Mariotte's law for perfect gases.
It came from observation that, keeping T constant, the measures (v, p) was closely aligned on an hyperbola branch, thus suggesting that p was of the form p=Constant/v.
But this relation is mostly written p*v = constant.

On the other side one keeps writitng P*v = n*R*T and not (as I said in my answer f(p, v, T) = 0) p*V-n*R*T = 0.

IMO, something like f(p, v, T) is less confusing and more general when we cannot isolate one thermodynamical quantity in terms of the two others.

Here is a (maybe) more formal way to do what you askd for:


# pressure p, volume v and temperature T are linked by an equation of state
# which, for Van der Walls gases takes ths form

f := (p, v, T) -> p - R*T/(v-b) - a/v^2

proc (p, v, T) options operator, arrow; p-R*T/(v-b)-a/v^2 end proc


# assuming volume v > covolume b > 0 one gets

normal(f(p, v, T) = 0);

numer(lhs(%)) = rhs(%)

(R*T*v^2+b*p*v^2-p*v^3-a*b+a*v)/((-v+b)*v^2) = 0


R*T*v^2+b*p*v^2-p*v^3-a*b+a*v = 0


# get the coefficients of v^n, n=0..3

CV := [seq(C[k], k in [seq](3..0, -1))] =~ [coeffs(lhs(%), v)]

[C[3] = -p, C[2] = R*T+b*p, C[1] = a, C[0] = -a*b]


# symbolic form

('f(p, v, T) = 0') &implies (add('C'[k]*v^k, k=0..3)=0)

`&implies`(f(p, v, T) = 0, v^3*C[3]+v^2*C[2]+v*C[1]+C[0] = 0)


# extended form

('f(p, v, T) = 0') &implies (eval(op(2, %), CV))

`&implies`(f(p, v, T) = 0, -p*v^3+(R*T+b*p)*v^2+a*v-a*b = 0)






To complete @Rouben Rostamian  advice:

sys := eval({eqn1, eqn2, eqn3, eqn4, i(0) = 11437, r(0) = 1077, s(0) = 1770000, t(0) = 1087}), {i(t), r(t), s(t), t(t)}:


will present each equation and initial condition on a different line to ease the visual analysis.

Also useful when womething goes wrong, is to check if the functions which appear in "diff" are also the functions which appear in initial or boundary conditions:

# functions you take derivatives
dfn := select(has, indets(sys, function), diff):
dfn := map2(op, 0, map2(op, 1, dfn));

# functions you define initial conditions for
fn0 := map2(op, 0, lhs~(ic));

Finally, check for all possibly non numerical parameters.
In your initial worksheet:

sys := {eqn1, eqn2, eqn3, eqn4, i(0) = 11437, r(0) = 1077, s(0) = 1770000, t(0) = 1087}:

# independent variable (t)
indepvar := op(2, select(has, indets(sys, function), diff)[1]);

# remaining parameters
indets(sys, name) minus {indepvar}
                            {T, mu}

Here you have two formal parameters (using print~(sys) reveals mu has not been set to 0.04 as intended).
As @Rouben Rostamian said it is unsafe to use mu and mu[1]; look at this

mu    := 2;
mu[1] := 3;

mu, mu[1];
                             mu, 3
mu[1] := 3;
mu := 2;
mu, mu[1];
                            2, 2[1]

Two possibilities to avoid this issue:

mu[0] := 2:
mu[1] := 3:
'mu[0]' = mu[0], 'mu[1]' = mu[1];

                      mu[0] = 2, mu[1] = 3

mu := 2:
mu__1:= 3:
'mu' = mu, 'mu__1' = mu__1;

                       mu = 2, mu__1 = 3

A last point, in case you really don't want to set T to some numerical value you can do this

dsn := dsolve(sys, numeric, parameters=[T]):

# Specialize dsn for some value of T

# Plot as usual
plots:-odeplot(dsn, ....)

One advantage of keepping free some parameters can be seen here:


dsn := dsolve({diff(x(t), t)=a, x(0)=0}, numeric, parameters=[a]):

# Compare the solution for different values of T
p := NULL:
for tau from 1 to 3 do
  p := p, plots:-odeplot(dsn, [t, x(t)], t=0..1, linestyle=tau)
end do:



one of my rare moments of inspiration

There are already many questions and answers about the SIR model.
Maybe you could start by looking at them, maybe you will find something interesting

In any case, it would be better if you would post your worksheet using the big green up arrow in the menu bar.

See you soon

2 3 4 5 6 7 8 Last Page 4 of 99