8 years, 67 days

## @zaidoun You're right, solz (ou...

@zaidoun

You're right, solz (output (5)) is the solution.

Please: note that the solution I obtained here was already found by Kitonum.

## @zaidoun The command substitut...

@zaidoun

The command

substitutes any occurences of z by exp(u), thus u becomes the "new variable".
The idea, following the initial remark Carl made, is to reduce the search range which extended to 5*10^10 to somethong smaller (z = 5*10^10 means u = log(5*10^10) ~25).
One problem is that your initial range for the imaginary part of z is -1000.. 5*10^10.
To cover this range with the change of variable z --> exp(u) is impossible. So, what I suggest you is:

1. solve F(z) for a reasonable range: previous answers showed RootFinding works well for re=-1000..1000 and im=-1000..5000.
2. Next solve G(u) on the extentendedrange  re=-1000..1000, im=log(5000)..log(5*10^10)

The command

enables recovering z solutions solz from u solutions solu.
Remain solu is a list, so the command exp(solu) won,'t work. To recover solu you need:

1. to build a loop over all the elements of solu and compute the exponential of each of them,
2. or to apply exp to all the elements of solin a single shot:
1. either by using the tilda symbol (see the corresponding help page for the tilda symbol)
2. or by writing solz := map(exp, solu)

Finally the two last command evaluate G(u) for each u in solu and F(z) for each z in solz.
And YES, the result is (then 4 times the same, a thing that ought to be verified)

## @zaidoun  Maybe you could find som...

@zaidoun

Maybe you could find some ideas here

 >
 >
 (1)
 >
 (2)
 >
 (3)
 (4)
 >
 (5)
 > eval~(G, u=~sol__u); eval~(F, z=~sol__z);
 (6)
 >

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