620 Reputation

2 Badges

3 years, 19 days

MaplePrimes Activity

These are replies submitted by mmcdara

@Rouben Rostamian  

Given that  Xbar = Sum(x[i], i=1..n)/n one easily finds that diff(Xbar, x[i]) = 1/n whatever the value of i
Then diff(b, x[i]) = diff(b, Xbar)*diff(Xbar, x[i]) = diff(b, Xbar)/n
(which is just "good sense" if we refer to the statistical meaning of Xbar).

Your procedure c_2 computes diff(b, x[i]). 
You will find that  sum(c_2(my_x, my_xbar, i), i=1..numelems(my_x)) = 0  ***(see below)
But, given diff(b, x[i]) = diff(b, Xbar)/n, applying the summation on the lhs and rhs of this equality should imply diff(b, Xbar)=0
You find a value (c_1) of -1.479...

This comes from the fact that, contrary to the requirement, x[i] and Xbar are not independent.
This hypothesis farazhedayati asked you to use is inconsistent

About the result sum(c_2(my_x, my_xbar, i), i=1..numelems(my_x)) = 0 
This result holds for any list my_x
(try for instance my_x := convert(Statistics:-Sample(Uniform(0, 1), 17), list) ).
This is a direct consequence of the fact that Xbar = Sum(x[i], i=1..n)/n


"Please consider Xbar and x[i] as independent variables" : Wrong, as xbar=sum(x[i], i=1..n)/n they are obviously not independent
diff( sum(x[i], i=1..n)/n, x[i]) = 1/n whatever the value of i


Sorry for this late reply.

Here is an attempt to answer your question.
Be carefull: the method I used in my previous post can give false results if the integrand is singukar over the integration domain (look yellow highlighted red text)



# some examples

  a := diff(f(x,y), x$3) + diff(f(x,y), [x, y]) + diff(f(x,y), y$2) + f(x, y) + diff(g(x, y), y) + h(y);
# a := diff(f(x,y), x$3) + diff(f(x,y), [x, y]);
# a := f(x, y) + h(y);

diff(diff(diff(f(x, y), x), x), x)+diff(diff(f(x, y), x), y)+diff(diff(f(x, y), y), y)+f(x, y)+diff(g(x, y), y)+h(y)


# define the alias ExI


# hasfun(a, diff) returns true if expression a contains terms with derivatives

if hasfun(a, diff) then
  b  := [op(select(has, a, diff))];                 # extract all terms with derivatives

  c  := op~(-1, [op(%)]);                           # extract the last differentiation variable for each term in b
                                                    # c will be used to integrate with respect to these last
                                                    # differentiation variables
  d  := map(C -> if C=x then y else x end if, c);   # d denotes the other integration variable

  zip((u, v) -> int(u, v), b, c);                   # integrate ech term in b along its last differentiation variable
  I1 := add(zip((u, v) -> int(u, v), %, d));        # integrate each previous term along the second integration variable
end if:

# add to I1 the double integrations with respect to x and y of all the terms that do not contain derivatives

I2 := I1 + ExI(ExI(select(`not`(has), a, diff), x), y)



int(diff(diff(f(x, y), x), x), y)+f(x, y)+int(diff(f(x, y), y), x)+int(g(x, y), x)+int(int(f(x, y), x), y)+x*(int(h(y), y))


# boundaries of the integration domain Omega (cartesian prodict of BX by BY)

BX := [x__1, x__2];
BY := [y__1, y__2];

[x__1, x__2]


[y__1, y__2]


# example

I3 := eval(I2, [f(x, y)=x*cos(y), g(x, y)=exp(x*y), h(y)=1/y^3])



# matrix of the evaluations of I3 at the 4 corners of Omega

M := Matrix(2, 2, (i, j) -> eval(I3, [x=BX[i], y=BY[j]]))

M := Matrix(2, 2, {(1, 1) = `#msub(mi("x"),mi("1"))`*cos(`#msub(mi("y"),mi("1"))`)+exp(`#msub(mi("x"),mi("1"))`*`#msub(mi("y"),mi("1"))`)/`#msub(mi("y"),mi("1"))`-(1/2)*`#msub(mi("x"),mi("1"))`/`#msub(mi("y"),mi("1"))`^2, (1, 2) = `#msub(mi("x"),mi("1"))`*cos(`#msub(mi("y"),mi("2"))`)+exp(`#msub(mi("x"),mi("1"))`*`#msub(mi("y"),mi("2"))`)/`#msub(mi("y"),mi("2"))`-(1/2)*`#msub(mi("x"),mi("1"))`/`#msub(mi("y"),mi("2"))`^2, (2, 1) = `#msub(mi("x"),mi("2"))`*cos(`#msub(mi("y"),mi("1"))`)+exp(`#msub(mi("x"),mi("2"))`*`#msub(mi("y"),mi("1"))`)/`#msub(mi("y"),mi("1"))`-(1/2)*`#msub(mi("x"),mi("2"))`/`#msub(mi("y"),mi("1"))`^2, (2, 2) = `#msub(mi("x"),mi("2"))`*cos(`#msub(mi("y"),mi("2"))`)+exp(`#msub(mi("x"),mi("2"))`*`#msub(mi("y"),mi("2"))`)/`#msub(mi("y"),mi("2"))`-(1/2)*`#msub(mi("x"),mi("2"))`/`#msub(mi("y"),mi("2"))`^2})


# definite integration of I3 over Omega

I4 := M[2, 2] + M[1, 1] - M[1, 2] - M[2, 1]



# verification

A := eval(a, [f(x, y)=x*cos(y), g(x, y)=exp(x*y), h(y)=1/y^3]);

# WATCHOUT : the above strategy doesn't handle possible singularities inside Omega
ia := int(int(A, x=BX[1]..BX[2]), y=BY[1]..BY[2]);

# I cheat here

IA := int(int(A, x=BX[1]..BX[2]), y=BY[1]..BY[2]) assuming BX[1] < BX[2], 0 < BY[1], BY[1] < BY[2]



Warning, unable to determine if 0 is between y__1 and y__2; try to use assumptions or use the AllSolutions option


int((x__1*sin(y)*y^3-x__2*sin(y)*y^3-exp(x__1*y)*x__1*y^2+exp(x__2*y)*x__2*y^2+exp(x__1*y)*y-exp(x__2*y)*y-x__1+x__2)/y^3, y = y__1 .. y__2)




simplify(IA - I4)




int(diff(1/x, x), x=-1..1);

# but

int(diff(F(x), x), x);
eval(%, F(x)=1/x);
eval(%, x=1) - eval(%, x=-1)














You're right


L := (r, theta) -> -2* diff(f(r, theta), r$2)+12*f(r, theta)-12*diff(f(r, theta), theta$2)/r+0.6*(diff(f(r, theta), theta, theta, r, r))*(1/4)+0.5*(diff(f(r, theta), theta, theta, theta, theta))*(1/4)

proc (r, theta) options operator, arrow; -2*(diff(f(r, theta), `$`(r, 2)))+12*f(r, theta)-12*(diff(f(r, theta), `$`(theta, 2)))/r+.6*(diff(f(r, theta), theta, theta, r, r))*(1/4)+.5*(diff(f(r, theta), theta, theta, theta, theta))*(1/4) end proc


S := L(r, theta);

-2*(diff(diff(f(r, theta), r), r))+12*f(r, theta)-12*(diff(diff(f(r, theta), theta), theta))/r+.1500000000*(diff(diff(diff(diff(f(r, theta), r), r), theta), theta))+.1250000000*(diff(diff(diff(diff(f(r, theta), theta), theta), theta), theta))


IS := int(int(op(1, S), r), theta) + int(int(op(2, S), r), theta) + int(int(op(3, S), theta), r) + int(int(op(4, S), theta), r) + int(int(op(5, S), theta), r)

int(-2*(diff(f(r, theta), r)), theta)+int(int(12*f(r, theta), r), theta)+int(-12*(diff(f(r, theta), theta))/r, r)+.1500000000*(diff(diff(f(r, theta), r), theta))+int(.1250000000*(diff(diff(diff(f(r, theta), theta), theta), theta)), r)


F := (sin(-4.700000000*10^(-6)+4.700000000*r)-0.1369508410*sinh(-4.700000000*10^(-6)+4.700000000*r))*cos(6*theta);



subs(f(r, theta)=F, IS):

# Exact primitive

ISF := eval(%);

.1666666667*(-9.400000000*cos(-0.4700000000e-5+4.700000000*r)+1.287337905*cosh(-0.4700000000e-5+4.700000000*r))*sin(6.*theta)+.1666666667*(-2.553191489*cos(-0.4700000000e-5+4.700000000*r)-.3496617217*cosh(-0.4700000000e-5+4.700000000*r))*sin(6.*theta)+(0.3384000000e-3+0.*I)*sin(6.*theta)*Ei(1., -(4.700000000*I)*r)+4.930207104*sin(6.*theta)*Ei(1., -4.700000000*r)-4.930253448*sin(6.*theta)*Ei(1., 4.700000000*r)+(-113.0973355-0.5315574770e-3*I)*sin(6.*theta)*csgn(r)+(72.00000000+0.3384000000e-3*I)*sin(6.*theta)*Si(4.700000000*r)-.9000000000*(4.700000000*cos(-0.4700000000e-5+4.700000000*r)-.6436689527*cosh(-0.4700000000e-5+4.700000000*r))*sin(6*theta)+0.2700000000e-7*sin(6.*theta)*(-212765957.4*cos(-0.4700000000e-5+4.700000000*r)-29138476.81*cosh(-0.4700000000e-5+4.700000000*r))


Ar := [0, 1]:
At := [0, 2*Pi]:

sISF := simplify(ISF):
M := Matrix(2, 2, (i,j) -> evalf(eval(limit(sISF, r=Ar[i]), theta=At[j])));

M := Matrix(2, 2, {(1, 1) = Float(undefined), (1, 2) = Float(undefined), (2, 1) = 0.+0.*I, (2, 2) = 0.+0.*I})


IL := M[1, 1] + M[2, 2] - M[1, 2] - M[2, 1]






Thanks, I will remember that


Probably not the simplest way to do this, but it works...

expr := 3*L*r/m/l/s+L/m:
MyRule := beta=L/m;

subs(L=solve(MyRule, L), expr);

eval(L[1, 1](r, theta), u__0(r, theta) = f[1, 1](r, theta));

# or

subs(u__0(r, theta) = f[1, 1](r, theta), L[1, 1](r, theta))

@Thomas Richard 

Answered from home, 

Thank you Thomas.

Since your optimization call succeeds both with method=ego and method=diffevol, there is probably no need to try other methods and options - apart from curiosity, which is okay. ;-)
Right, but the curiosity wasn't accidental: In my real application the solution of the ODE system depends on 10 parameters
(with quite different ranges of variation, but I normalized both of them between 0 and 1 to avoid numerical difficulties).
The first time I used GOT with diffevol I got an error message concerning the time limit.  Not the message you find in GetLastSolution help page, but an error message (I will be able to send it tomorrow) wich, from memory, said "error in convert/radical ... in < here some procedure, not always the same > ".
So I began to try other things.
Finally the problem was solved by setting the option "localrefinement" (sorry, I cite from memory) to false and by using diffevol ("ego" is much slower, so the idea to reduce the size of the design of experiments to build the response surfaces)

If you are interested I could post the true application tomorrow.


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

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


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)


Maybe you could find some ideas here


restart; Digits := 20

r1 := sqrt(2^2-I*z);



r2 := sqrt(2^2+I*z);



F := -r2*(r1^2-(1/5)*2^2)^2*cosh((1/150)*r1*Pi)*sinh((1/150)*r2*Pi)+r1*(r2^2-(1/5)*2^2)^2*cosh((1/150)*r2*Pi)*sinh((1/150)*r1*Pi)+z*(-r1^2+r2^2)*cosh((1/150)*r1*Pi)*cosh((1/150)*r2*Pi)





G := subs(z = exp(u), F);



[3.8593940676611702887+21.991148575128552669*I, 3.8593940676611702887+15.707963267948966192*I, 3.8593940676611702888+9.4247779607693797150*I, 3.8593940676611702887+3.1415926535897932385*I]


[-47.436599289272340540+0.11313803637424101255e-16*I, -47.436599289272340540+0.14857945353770406688e-16*I, -47.436599289272340545+0.18402087070116712123e-16*I, -47.436599289272340540-0.17720708581731527164e-17*I]


eval~(G, u=~sol__u);
eval~(F, z=~sol__z);

[-0.10688269926925727398e-14+0.1e-15*I, -0.14361096901711096670e-14+0.*I, -0.17178880581261977346e-14-0.7e-15*I, 0.10314134873926846362e-15+0.1e-15*I]


[-0.10688269926925727398e-14+0.1e-15*I, -0.14361096901711096670e-14+0.*I, -0.17178880581261977346e-14-0.7e-15*I, 0.10314134873926846362e-15+0.1e-15*I]





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



sys := [diff(Y(x, t), t, t) = 0, Y(x, 0) = 0, (D[2](Y))(x, 1) = 0];

`Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895`


[diff(diff(Y(x, t), t), t) = 0, Y(x, 0) = 0, (D[2](Y))(x, 1) = 0]


Y(x, t) = 0





@Carl Love 

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

  • I use to visit the site 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



    • 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, these alse prompts "are a total nuisance to remove"


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.

1 2 3 4 5 6 7 Last Page 1 of 28