8 years, 52 days

Don't have Maple 2018 right now,&nbs...

Don't have Maple 2018 right now, here is the result obtained with Maple 2015.2

 > restart:
 > interface(version); sys := [diff(Y(x, t), t, t) = 0, Y(x, 0) = 0, (D[2](Y))(x, 1) = 0]; pdsolve(sys);
 (1)
 >

@Carl Love  It's the kind of t...

It's the kind of thing you can do when there's nothing more interesting around here

• I use to visit the site https://oeis.org/?language=english and I went to search if the sequence
1, 28, 110, 384, 1316, 4496, 15352, 52416, 178960, 611008 ...  would not have interesting properties.
Unfortunately this sequence is not identified, but, if you set f[1]=1 and f[2]=M, then f[n] = A[n]*M+B[n] where
A = 0, 1, 4, 14, 48, 164, 560 ... and B = 1, 0, -2, -8, -28, -96, -328 ...
Sequence A[n], from n> 1, is dubbed A007070
Sequence B[n], from n> 2, is dubbed A106731
They are both related and many characterizations for both of them, including generating functions and expansions of rational fractions are given.
Unfortunately I wasn't capable to find any clue to prove the sequence generated with M=28 doesn't contain any square (except f[1] of course).
Maybe you would be more successful?

• So I had a little fun, here are a few results I got

Observations:

• For f[1]=1 only 25 values of f[2] in the range [2..675] contains (at least) one f[n] which is a square when n varies from 1 to 1000

• The corresponding value of n is rather small... maybe it is even unique?

• One also generate sequences with f[1]=p and f[2]=m (m>p). One finds that for some couples (m, p) there exist different values of n such that f[n] is a square

@Carl Love  Thanks Carl. I agree, ...

@Carl Love

Thanks Carl.
I agree, these alse prompts "are a total nuisance to remove"

@phiost No problem.In case you woul...

@phiost

No problem.

In case you would be nevertheless interested in coloring "3D level curves" and identifying them by a legend you could find a few ideas here.

contourplot2.mw

@Carl Love  Just a remark about th...

@Carl Love

Just a remark about the "prompt",

If you copy-paste your first conditional sequence into a worksheet (worksheet mode) then its five line are preceeded with a prompt.
But, as they are in the same single execution group, the result is the one expected.
The prompts obtained by this copy-paste "are the same" as those you get when applying menu>edit>join on multiple groups.
Then it seems there are two kind of prompts, saying the "active" ones which begin a group and the "passive" ones inside a group.
So maybe it's notcompletely true to say "in other words, a single command prompt".

@vv  Smart, as Iwrote to silvergas...

@vv

Smart, as Iwrote to silvergasshoper I had done that years ago.
Your solution is definitely better, I thumb up.

A point of detail...In the case where th...

A point of detail...

In the case where the thermal conductivity is discontinuous, your scheme doesn't guarantee the (physical) continuity of the thermal flux (unless the mesh size is adpated in a correct way left and right of the discontinuity points).
The classical method to avoid such unphysical solutions consists in deriving finite difference approximations that explicitely account for the continuity equations for the temperature and the thermal flux.
One can show that the correct thermal continuity at its points of disconinuity is the harmonic mean of its values left and right to these points.

On coarse grids the "harmonic" method is by far the better. On finer grids its superiority on the classical method (the one you used) tends to lessen. More of this the "harmonic" method is more time consuming and may be found complex to implement efficiently on unstructured meshes.

 > restart:
 > phi := x -> lambda(x)*diff(theta(x), x)
 (1)
 > # let x = a a discontinuity point for lambda, # for instance lambda := x -> piecewise(x< a, lambda__L, lambda__R)
 (2)
 > # for the heat equation, continuity equations for # theta and flux phi  must be written : theta__continuity := limit(theta(x), x=a, left) = limit(theta(x), x=a, right); phi__continuity   := limit(phi(x), x=a, left) = limit(phi(x), x=a, right);
 (3)
 > # finite difference approximation approximation : # # Let x__1 < .. < x__(i-1) < a=x__i < x__(i+1) < .. < x__N a cpllection of # x points. # Let omega__(k+1/2) the finite volume [x__k, x__(k+1)]. # # In finite volume approaches a staggered mesh is used where temperatures # are estimated at the mid points x__(k+1/2) (k=1..N-1) and fluxes phi are # estimated at the finite volumes boundaries x__k (k=1..N). # # Let DL the divided difference the approximation of phi(x) in omega__(i-1/2) # # Let p = i-1/2 and x__i__L = limit(x__i - epsilon, epsilon=0^+) DL := lambda__L * ( theta(x__i__L) - theta(x__p) ) / ( x__i__L - x__p )
 (4)
 > # Identically let DR the divided difference the approximation of phi(x) in omega__(i+1/2) # # Let q = i+1/2 and x__i__R = limit(x__i + epsilon, epsilon=0^+) DR := lambda__R * ( theta(x__q) - theta(x__i__R) ) / ( x__q - x__i__R )
 (5)
 > # Remind theta is defined as points x__p and x__q alone. # This means x__i is undefined and has to be estimated. # This approximation relue supon the two continuity equations : # by continuity of theta at point x__i one has cond1L := theta(x__i__L) = theta(x__i); cond1R := theta(x__i__R) = theta(x__i); # by continuity of phi at point x__i one has cond2 := DL = DR; # at last, obviously, cond3 := x__i__L = x__i; cond4 := x__i__R = x__i;
 (6)
 > COND := subs({cond1L, cond1R, cond3, cond4}, cond2)
 (7)
 > # for the sake of simplicity, let's assume x__i-x__p = x__q-x__i = h/2 mesh := {x__p = x__i - h/2, x__q = x__i + h/2}
 (8)
 > solve(COND, theta(x__i))
 (9)
 > subs(mesh, %);
 (10)
 > # temperature at point a (or at any point x__k where k is now an integer) theta__a := simplify(%);
 (11)
 > # Injecting theta__a into DL (could be done with DR too) gives : subs(theta(x__i__L)=theta__a, DL); subs(cond3, %); subs(mesh, %); # thermal flux at point x=a (or at any point x__k where k is now an integer) # (also known as "interfacial thermal flux") phi__a := simplify(%);
 (12)
 > # The equivalent interfacial conductivity lambda__eq is de fined by the relation : rel := lambda__eq*( theta(x__i+(1/2)*h) - theta(x__i-(1/2)*h) ) / h = phi__a: Lambda__eq := solve(rel, lambda__eq);
 (13)
 > # if lambda(x) is continuous at point x = a one has subs(lambda__R = lambda__L, Lambda__eq);
 (14)
 > # When continuity doesn't hold the equivalent interfacial conductivity is twice the # harmonic mean of the left and right conductivities. # # In any case the correct divided differences approximation of the thermal flux is # the expression phi__a above: only this one guarantees the continuity of the temperature # and of the thermal flux at points where the thermal conductivity is discontinuous.
 >

To complete Tomleslie's answer, you ...

To complete Tomleslie's answer, you could be interested to know that the value of  time "t0" can also be found by using the events feature of dsolve
(type "events" in the help and look to the  dsolve, numeric, Events (events) help page)

With tomleslies' notations ("events" are supported by rkf45 [default}, ck45 and rosenbrock methods)

sol1:= dsolve(
{ sys, x(0)=0, y(0)=0, z(0)=cos(1), w(0)=sin(1)},
type=numeric, method=rosenbrock, events=[[[x(t)<0, And(-y(t)<-5.5, y(t) < 6.5)], halt]]
);

LargeEnoughTime := 20;  # set LargeEnoughTime to a value larger than the one returned by the warning message
sol1(LargeEnoughTime):

EndTime := op(sol1(eventfired=[1]));
sol1(EndTime);

EndTime := 7.84831340510192

@Carl Love  By the way, what is th...

@Carl Love

By the way, what is the risk to write simultaneously omega and omega[0]?
Can't that fool Maple?

Hi,Here is a lengthy method to obtain th...

Hi,
Here is a lengthy method to obtain the result by hand.
I don't say a more astute solution could not be found

 > restart:
 > expr_1 := 512*b^9 + (2303*a + 2304)*b^8 + (4616*a^2 + 9216*a + 4608)*b^7 + (5348*a^3 + 16128*a^2 + 16128*a + 5376)*b^6
 > + (4088*a^4 + 16128*a^3 + 24192*a^2 + 16128*a + 4032)*b^5 + (1946*a^5 + 10080*a^4 + 20160*a^3 + 20160*a^2
 > + 10080*a + 2016)*b^4 + (728*a^6 + 4032*a^5 + 10080*a^4 + 13440*a^3 + 10080*a^2 + 4032*a + 672)*b^3
 > + (116*a^7 + 1008*a^6 + 3024*a^5 + 5040*a^4 + 5040*a^3 + 3024*a^2 + 1008*a + 144)*b^2
 > + (26*a^8 + 144*a^7 + 504*a^6 + 1008*a^5 + 1260*a^4 + 1008*a^3 + 504*a^2 + 144*a + 18)*b + 9*a^8
 > + 36*a^7 + 84*a^6 + 126*a^5 + 126*a^4 + 84*a^3 + 36*a^2 + 9*a + 1
 >
 (1)
 > # to found : expr_0 := (a+2*b+1)^9 - a*(a-b)^8
 (2)
 > # just to check simplify(expand(expr_0)-expand(expr_1));
 (3)
 > # first visual observation c_0 := coeff(expr_1, b, 0); expand((a+1)^9); # thus # c_0 := (a+1)^9 - a^9 : # check expand(c_0 - ((a+1)^9 - a^9) ); c_0 := (a+1)^9 - a^9;
 (4)
 > # could it be that coeff(expr_1, b, 1) is of a comparable form? c_1  := coeff(expr_1, b, 1); test := expand((a+1)^8); for n from 0 to 8 do    coeff(c_1, a, n) / coeff(test, a, n) end do;
 (5)
 > # thus : # c_1 := 18*(a+1)^8 + (26-18)*a^8: # check expand(c_1 - (18*(a+1)^8 + (26-18)*a^8) ); c_1 := 18*(a+1)^8 + (26-18)*a^8;
 (6)
 > # coeff(expr_1, b, 2) c_2  := coeff(expr_1, b, 2); test := expand((a+1)^7); for n from 0 to 7 do    coeff(c_2, a, n) / coeff(test, a, n) end do;
 (7)
 > # thus : # c_2 := 144*(a+1)^7 + (116-144)*a^7: # check expand(c_2 - (144*(a+1)^7 + (116-144)*a^7) ); c_2 := 144*(a+1)^7 + (116-144)*a^7;
 (8)
 > # coeff(expr_1, b, 3) c_3  := coeff(expr_1, b, 3); test := expand((a+1)^6); for n from 0 to 6 do    coeff(c_3, a, n) / coeff(test, a, n) end do; # thus : # c_3 := 672*(a+1)^6 + (728-672)*a^6: # check expand(c_3 - (672*(a+1)^6 + (728-672)*a^6) ); c_3 := 672*(a+1)^6 + (728-672)*a^6;
 (9)
 > # coeff(expr_1, b, 4) c_4  := coeff(expr_1, b, 4); test := expand((a+1)^5); for n from 0 to 5 do    coeff(c_4, a, n) / coeff(test, a, n) end do; # thus : # c_4 := 2016*(a+1)^5 + (1946-2016)*a^5: # check expand(c_4 - (2016*(a+1)^5 + (1946-2016)*a^5) ); c_4 := 2016*(a+1)^5 + (1946-2016)*a^5;
 (10)
 > # coeff(expr_1, b, 5) c_5  := coeff(expr_1, b, 5); test := expand((a+1)^4); for n from 0 to 4 do    coeff(c_5, a, n) / coeff(test, a, n) end do; # thus : # c_5 := 4032*(a+1)^4 + (4088-4032)*a^4: # check expand(c_5 - (4032*(a+1)^4 + (4088-4032)*a^4) ); c_5 := 4032*(a+1)^4 + (4088-4032)*a^4;
 (11)
 > # coeff(expr_1, b, 6) c_6  := coeff(expr_1, b, 6); test := expand((a+1)^3); for n from 0 to 3 do    coeff(c_6, a, n) / coeff(test, a, n) end do; # thus : # c_6 := 5376*(a+1)^3 + (5348-5376)*a^3: # check expand(c_6 - (5376*(a+1)^3 + (5348-5376)*a^3) ); c_6 := 5376*(a+1)^3 + (5348-5376)*a^3;
 (12)
 > # coeff(expr_1, b, 7) c_7  := coeff(expr_1, b, 7); test := expand((a+1)^2); for n from 0 to 2 do    coeff(c_7, a, n) / coeff(test, a, n) end do; # thus : # c_7 := 4608*(a+1)^2 + (4616-4608)*a^2: # check expand(c_7 - (4608*(a+1)^2 + (4616-4608)*a^2) ); c_7 := 4608*(a+1)^2 + (4616-4608)*a^2;
 (13)
 > # coeff(expr_1, b, 8) c_8  := coeff(expr_1, b, 8); test := expand((a+1)); for n from 0 to 1 do    coeff(c_8, a, n) / coeff(test, a, n) end do; # thus : # c_8 := 2304*(a+1) + (2303-2304)*a: # check expand(c_8 - (2304*(a+1) + (2303-2304)*a) ); c_8 := 2304*(a+1) + (2303-2304)*a;
 (14)
 > # résumé # Remark c_8 is of the form 2304*(a-1) - 1 seq(print(c_||k), k=0..9);
 (15)
 > # Observe the coefficietnts of the second term : C := [ -1, 8, -28, 56, -70, 56, -28, 8, -1 ]
 (16)
 > # these are the coefficients of the powers of a in -(a-1)^8: expand(-(a-1)^8); seq(coeff(expand(-(a-1)^8), a, k), k=0..8)
 (17)
 > # Thus the sequence of all the op(2, c_k) are the terms of expand(-a*(a-1)^8);
 (18)
 > # Then expr_0 is of the form # # 512*b^9 + add(c_||k*b^k, k=0..8) = 512*b^9 + add(op(1, c_||k)*b^k, k=0..8) + add(op(2, c_||k)*b^k, k=0..8) #                                  = 512*b^9 + add(op(1, c_||k)*b^k, k=0..8) + add(coeff(T, a, 9-k)a^(9-k)*b^k, k=0..8) #                                  = 512*b^9 + add(op(1, c_||k)*b^k, k=0..8) - a*(a-b)^8 T1 := - a*(a-b)^8
 (19)
 > T2 := simplify( expr_0 - T1);
 (20)
 > simplify(expr_0 - (T1+T2))
 (21)
 >

@chandrashekhar  I did again a cop...

@chandrashekhar

I did again a copy-paste of your initial expression and of its manual simplification of yours.
I used, as you suggested, simplify(s1-yourexpression,size): this doesn't give 0.
Two possibilities:

1. the copy paste did not work well (but I have serious doubt about that)
2. you played me when writting that simplification(s1-yourexpression,size) gives 0!
First of all, simplification IS NOT A MAPLE COMMAND ... but maybe it was a typo ???

I really not like this kind of dishonesty.

 > restart:
 > S0 := 4*beta*c*lambda^2*mu^2+2*beta*c*lambda*mu^3+2*beta*c*lambda*mu^2*nu+beta*c*mu^3*nu-4*c*lambda^2*mu^2*sigma-2*c*lambda*mu^3*sigma-2*c*lambda*mu^2*nu*sigma-c*mu^3*nu*sigma+4*beta*lambda^3*sigma+4*beta*lambda^2*mu*sigma+2*beta*lambda^2*nu*sigma+2*beta*lambda*mu^2*sigma+2*beta*lambda*mu*nu*sigma+beta*mu^3*sigma+beta*mu^2*nu*sigma+4*lambda^2*mu^2*sigma+2*lambda*mu^3*sigma+2*lambda*mu^2*nu*sigma+mu^3*nu*sigma
 (1)
 > S1 := simplify(S0, size);
 (2)
 > YourExpression := beta*mu^3+(2*lambda+nu)*(beta*mu^2+2*beta*lambda*mu+2*beta*lambda^2+(1-c)*(2*lambda+mu)*mu^2)
 (3)
 > simplify(S1-YourExpression, size)
 (4)
 > # simplification IS NORT A MAPLE COMMAND simplification(S1-YourExpression, size)
 (5)
 > vars := indets(S0, name);
 (6)
 > M := Matrix(numelems(vars), 3): k := 0: for v in vars do    k       := k+1:    M[k, 1] := degree(collect(S0, v), v);    M[k, 2] := degree(collect(S1, v), v);    M[k, 3] := degree(collect(YourExpression, v), v); end do: "Higher degree in each indet for S0, S1 and your manual simplification"; <   Vector[column](numelems(vars)+1, ["indet", vars[]])   |   <     Vector[row](3, ["Initial expression", "Maple simplified", "Manually simplified"]),     M   > >;
 (7)
 >
 >

@acer  I just copied-pasted myself...

@acer

I just copied-pasted myself (maybe I missed something to) but I found that the Maple simplification of the original expression XWAS NOT equivalent to the manually simplified one:

 > restart:
 > S0 := 4*beta*c*lambda^2*mu^2+2*beta*c*lambda*mu^3+2*beta*c*lambda*mu^2*nu+beta*c*mu^3*nu-4*c*lambda^2*mu^2*sigma-2*c*lambda*mu^3*sigma-2*c*lambda*mu^2*nu*sigma-c*mu^3*nu*sigma+4*beta*lambda^3*sigma+4*beta*lambda^2*mu*sigma+2*beta*lambda^2*nu*sigma+2*beta*lambda*mu^2*sigma+2*beta*lambda*mu*nu*sigma+beta*mu^3*sigma+beta*mu^2*nu*sigma+4*lambda^2*mu^2*sigma+2*lambda*mu^3*sigma+2*lambda*mu^2*nu*sigma+mu^3*nu*sigma
 (1)
 > S1 := simplify(S0, size);
 (2)
 > YourExpression := beta*mu^3+(2*lambda+nu)*(beta*mu^2+2*beta*lambda*mu+2*beta*lambda^2+(1-c)*(2*lambda+mu)*mu^2)
 (3)
 > simplify(S1-YourExpression)
 (4)
 >

@tomleslie  Thanks fot the check. ...

@tomleslie

Thanks fot the check.
I forget this problem and consider it as closed.

@tomleslie  I sincerely thank you ...

@tomleslie

I sincerely thank you for this very thorough analysis, which probably took you some time.
Does that enlighten me? I don't know, but it makes me feel a little more comfortable because I was sincerely beginning to wonder if I hadn't done something so stupid and I couldn't see it.

I guess if no one on the development team takes care of this problem, it will remain a quirk.

Just one question, being currently on holiday, I can only use Maple 2015: do you observe the same phenomenon with a more recent version, for example Maple 2018 or 2019?
If this were the case, I think it would be worthwhile to open a new thread on this "toto"; otherwise, we could consider that bygones are bygones

And last, "The anomaly is much too subtle for me to work out WTF is going on - sorry:-(", don't be sorry, you really did a fine job

@tomleslie

After verifications I have no such .bak files in the directory where I try to save the worksheet.
Going deeply  I observed that

• dodgy0 can be saved without any problems after having been executed
dodgy0.mw
Curious thing : there is no ouput returned after the 3 dsolve commands... there is something here that doesn't sound normal

• dodgy1 can be saved without any problems ONLY if it has not been executed (!!! icon)
dodgy1.mw
If you load  and execute it, you will see that each dsolve command nof return the classical proc(...) ...end proc output
But once executed I'm not able to save it (cmd+s, on Mac OS)

• dodgy2 can be saved without any problems after having been executed
dodgy2.mw
Its only difference with dodgy1: dsolve(..., method=rosenbrock) has been removed

• dodgy3 can be saved without any problems ONLY if it has not been executed (!!! icon)
dodgy3.mw
Clearly dsolve(..., method=rosenbrock), the only dsolve present in this worksheet, seems to be the cause.

PS: if the ODE has real coefficients instead of complex ones, there is no problem at all

If you can't reproduce those behaviours, maybe they come from some conflict between Maple 2015.2 and  Mac OS Mojave 10.14.3 ?

 First 105 106 107 108 109 110 111 Last Page 107 of 134
﻿