dharr

Dr. David Harrington

8270 Reputation

22 Badges

20 years, 359 days
University of Victoria
Professor or university staff
Victoria, British Columbia, Canada

Social Networks and Content at Maplesoft.com

Maple Application Center
I am a retired professor of chemistry at the University of Victoria, BC, Canada. My research areas are electrochemistry and surface science. I have been a user of Maple since about 1990.

MaplePrimes Activity


These are answers submitted by dharr

Inverse Laplace transforms often require assumptions to work. Here assuming alpha is positive works, at least for Ex 1.

Ex1.mw

file menu, bottom right corner - labeled options.

@salim-barzani I didn't use your ugly code because you did not have the correct x,t dependence in it at the time, though you since fixed that. But I also wanted general code - there is no reason to have special P,Q,R for every case when you can go directly to A[n].

Thanks for pointing out the missing factorial. You may be suggesting something else is missing but if so I didn't understand what you meant. So I think this code is correct now. The conclusion is now the same as @acer's in the thread here, i.e., that the paper has an error.

Laplace Adomian Decomposition method.

restart

with(Physics); with(inttrans); PDEtools:-declare(U(x, t), quiet)

pde0 := I*(diff(U(x, t), t))+diff(U(x, t), `$`(x, 2))+2*(diff(U(x, t)*conjugate(U(x, t)), x))*U(x, t)+U(x, t)^2*conjugate(U(x, t))^2*U(x, t)

I*(diff(U(x, t), t))+diff(diff(U(x, t), x), x)+(2*(diff(U(x, t), x))*conjugate(U(x, t))+2*U(x, t)*(diff(conjugate(U(x, t)), x)))*U(x, t)+U(x, t)^3*conjugate(U(x, t))^2

pde := expand(-I*pde0)

diff(U(x, t), t)-I*(diff(diff(U(x, t), x), x))-(2*I)*U(x, t)*(diff(U(x, t), x))*conjugate(U(x, t))-(2*I)*U(x, t)^2*(diff(conjugate(U(x, t)), x))-I*U(x, t)^3*conjugate(U(x, t))^2

Operator parts

LLu := diff(U(x, t), t); RRu := I*(diff(U(x, t), x, x)); NNu := -pde+LLu-RRu

diff(U(x, t), t)

I*(diff(diff(U(x, t), x), x))

(2*I)*U(x, t)*(diff(U(x, t), x))*conjugate(U(x, t))+(2*I)*U(x, t)^2*(diff(conjugate(U(x, t)), x))+I*U(x, t)^3*conjugate(U(x, t))^2

We treat each [U(x, t), diff(U(x, t), x), diff(conjugate(U(x, t)), x), conjugate(U(x, t))]like a simple variable v[1],..,v[4].

vars := [op(`minus`(indets(NNu), {t, x}))]; nvars := nops(vars); vvars := [seq(v[i], i = 1 .. nvars)]; NNv := subs(`~`[`=`](vars, vvars), NNu)

[U(x, t), diff(U(x, t), x), diff(conjugate(U(x, t)), x), conjugate(U(x, t))]

4

[v[1], v[2], v[3], v[4]]

(2*I)*v[1]*v[2]*v[4]+(2*I)*v[1]^2*v[3]+I*v[1]^3*v[4]^2

u[0] starts as the initial condition, form which we find A[0] and Ru[0]. vals[0] are the values of the 4 vvars corresponding to u[0]

u[0] := beta*exp(I*x); vals[0] := `assuming`([simplify(eval(subs(conjugate(U(x, t)) = conjugate(u[0]), U(x, t) = u[0], vars)))], [real]); Asym[0] := subs({seq(v[j] = vars[j][0], j = 1 .. nvars)}, NNv); A[0] := `assuming`([simplify(eval(NNv, `~`[`=`](vvars, vals[0])))], [real]); Ru[0] := `assuming`([simplify(eval(RRu, U(x, t) = u[0]))], [real])

beta*exp(I*x)

[beta*exp(I*x), I*beta*exp(I*x), -I*beta*exp(-I*x), beta*exp(-I*x)]

(2*I)*U(x, t)[0]*(diff(U(x, t), x))[0]*conjugate(U(x, t))[0]+(2*I)*U(x, t)[0]^2*(diff(conjugate(U(x, t)), x))[0]+I*U(x, t)[0]^3*conjugate(U(x, t))[0]^2

I*beta^5*exp(I*x)

-I*beta*exp(I*x)

Next iteration.
u[1] becomes a function of t for the next Laplace transform. Asym is redundant here, jus to check.  At least the last term here and the second term in the paper (only ones with u[0]^3) are different.
A[1] here does not agree with the result in the numerical example.

n := 1; u[n] := invlaplace(laplace(Ru[n-1]+A[n-1], t, s)/s, s, t); vals[n] := `assuming`([simplify(eval(subs(conjugate(U(x, t)) = conjugate(u[n]), U(x, t) = u[n], vars)))], [real]); Asym[n] := eval((diff(subs({seq(v[j] = add(vars[j][i]*lambda^i, i = 0 .. n), j = 1 .. nvars)}, NNv), `$`(lambda, n)))/factorial(n), lambda = 0); A[n] := simplify(eval((diff(subs({seq(v[j] = add(vals[i][j]*lambda^i, i = 0 .. n), j = 1 .. nvars)}, NNv), `$`(lambda, n)))/factorial(n), lambda = 0)); Ru[n] := `assuming`([simplify(eval(RRu, U(x, t) = u[n]))], [real])

1

I*exp(I*x)*(beta^5-beta)*t

[I*exp(I*x)*beta*(beta^4-1)*t, t*(-beta^5+beta)*exp(I*x), t*(-beta^5+beta)*exp(-I*x), -I*exp(-I*x)*beta*(beta^4-1)*t]

(2*I)*U(x, t)[1]*(diff(U(x, t), x))[0]*conjugate(U(x, t))[0]+(2*I)*U(x, t)[0]*(diff(U(x, t), x))[1]*conjugate(U(x, t))[0]+(2*I)*U(x, t)[0]*(diff(U(x, t), x))[0]*conjugate(U(x, t))[1]+(4*I)*U(x, t)[0]*(diff(conjugate(U(x, t)), x))[0]*U(x, t)[1]+(2*I)*U(x, t)[0]^2*(diff(conjugate(U(x, t)), x))[1]+(3*I)*U(x, t)[0]^2*conjugate(U(x, t))[0]^2*U(x, t)[1]+(2*I)*U(x, t)[0]^3*conjugate(U(x, t))[0]*conjugate(U(x, t))[1]

exp(I*x)*t*(-beta^9+beta^5)

exp(I*x)*beta*(beta^4-1)*t

In a loop it is quite fast

for n to 5 do u[n] := invlaplace(laplace(Ru[n-1]+A[n-1], t, s)/s, s, t); vals[n] := `assuming`([simplify(eval(subs(conjugate(U(x, t)) = conjugate(u[n]), U(x, t) = u[n], vars)))], [real]); A[n] := simplify(eval((diff(subs({seq(v[j] = add(vals[i][j]*lambda^i, i = 0 .. n), j = 1 .. nvars)}, NNv), `$`(lambda, n)))/factorial(n), lambda = 0)); Ru[n] := `assuming`([simplify(eval(RRu, U(x, t) = u[n]))], [real]) end do

I*exp(I*x)*(beta^5-beta)*t

[I*exp(I*x)*beta*(beta^4-1)*t, t*(-beta^5+beta)*exp(I*x), t*(-beta^5+beta)*exp(-I*x), -I*exp(-I*x)*beta*(beta^4-1)*t]

exp(I*x)*t*(-beta^9+beta^5)

exp(I*x)*beta*(beta^4-1)*t

-(1/2)*exp(I*x)*(beta^9-2*beta^5+beta)*t^2

[-(1/2)*exp(I*x)*beta*(beta^8-2*beta^4+1)*t^2, -((1/2)*I)*exp(I*x)*beta*(beta^8-2*beta^4+1)*t^2, ((1/2)*I)*exp(-I*x)*beta*(beta^8-2*beta^4+1)*t^2, -(1/2)*exp(-I*x)*beta*(beta^8-2*beta^4+1)*t^2]

-((1/2)*I)*exp(I*x)*t^2*beta^5*(beta^8-2*beta^4+1)

((1/2)*I)*exp(I*x)*beta*(beta^8-2*beta^4+1)*t^2

-((1/6)*I)*exp(I*x)*(beta^4-1)^3*beta*t^3

[-((1/6)*I)*exp(I*x)*(beta^4-1)^3*beta*t^3, (1/6)*exp(I*x)*(beta^4-1)^3*beta*t^3, (1/6)*exp(-I*x)*(beta^4-1)^3*beta*t^3, ((1/6)*I)*exp(-I*x)*(beta^4-1)^3*beta*t^3]

(1/6)*exp(I*x)*t^3*beta^5*(beta^12-3*beta^8+3*beta^4-1)

-(1/6)*exp(I*x)*(beta^4-1)^3*beta*t^3

(1/24)*exp(I*x)*(beta^4-1)^4*beta*t^4

[(1/24)*exp(I*x)*(beta^4-1)^4*beta*t^4, ((1/24)*I)*exp(I*x)*(beta^4-1)^4*beta*t^4, -((1/24)*I)*exp(-I*x)*(beta^4-1)^4*beta*t^4, (1/24)*exp(-I*x)*(beta^4-1)^4*beta*t^4]

((1/24)*I)*exp(I*x)*t^4*(beta^2+1)^4*(beta-1)^4*beta^5*(beta+1)^4

-((1/24)*I)*exp(I*x)*(beta^4-1)^4*beta*t^4

((1/120)*I)*exp(I*x)*(beta^4-1)^5*beta*t^5

[((1/120)*I)*exp(I*x)*(beta^4-1)^5*beta*t^5, -(1/120)*exp(I*x)*(beta^4-1)^5*beta*t^5, -(1/120)*exp(-I*x)*(beta^4-1)^5*beta*t^5, -((1/120)*I)*exp(-I*x)*(beta^4-1)^5*beta*t^5]

-(1/120)*exp(I*x)*t^5*beta^5*(beta-1)^5*(beta+1)^5*(beta^2+1)^5

(1/120)*exp(I*x)*(beta^4-1)^5*beta*t^5

 

NULL

Download Laplace_Adomiian_Decomposition_method_4.mw

The error message indicates a problem with floats. Ideally, you should do the whole thing symbolically and then substitute in values, but that takes too long, probably because the expressions get too complicated. So you can convert the float parameters to rational; run through the loop and then evalf the results. Then you will need to change the parameters and run it again. It's possible that with numerical parameters you don't need the assumptions and/or the simplify (one or both), you you could try that for faster operation.

error-problem2.mw

It means there are some assumptions involving these variables. You can avoid displaying them with Options->Display->Assumed variables and choose No Annotation (or you might prefer Phrase).

There is no BlockDiagonalMatrix command. When the output of a command has the same command in it, that is a hint that the command does not exist. However, DiagonalMatrix does what you want.

Y8 := DiagonalMatrix([Y4, Y4_R]);

I don't understand the logic of what you are trying to do well enough to answer your other questions.

"f(p,n)= "the number of partions of a positive integer n into p parts. possible parts are only 0,1,4,6,8,9. these parts can be used several times. the order is not important.

restart

with(combstruct)

A word w is a multiset of letters (can have repeats, order doesn't matter) from s0,s1,s4,s6,s8,s9.

sys := {letters = Union(s0, s1, s4, s6, s8, s9), s0 = Atom, s1 = Atom, s4 = Atom, s6 = Atom, s8 = Atom, s9 = Atom, w = Set(letters)}

{letters = Union(s0, s1, s4, s6, s8, s9), s0 = Atom, s1 = Atom, s4 = Atom, s6 = Atom, s8 = Atom, s9 = Atom, w = Set(letters)}

Attribute grammar to calculate the sum s of the letters/numbers in the word

addition := {s(s0) = 0, s(s1) = 1, s(s4) = 4, s(s6) = 6, s(s8) = 8, s(s9) = 9, s(w) = Set(s(letters))}

{s(s0) = 0, s(s1) = 1, s(s4) = 4, s(s6) = 6, s(s8) = 8, s(s9) = 9, s(w) = Set(s(letters))}

Let z mark the number of letters (parts) in a word, and let u mark the value (sum).
The first few terms in the multivariate generating function are:

NULL
Order := 4; ser := agfseries(sys, addition, unlabeled, z, [[u, s]])[w(z, u)]

series(1+(u^9+u^8+u^6+u^4+u+1)*z+(u^18+u^17+u^16+u^15+u^14+u^13+2*u^12+2*u^10+2*u^9+2*u^8+u^7+u^6+u^5+u^4+u^2+u+1)*z^2+(u^27+u^26+u^25+2*u^24+u^23+2*u^22+2*u^21+2*u^20+2*u^19+4*u^18+3*u^17+4*u^16+2*u^15+3*u^14+3*u^13+3*u^12+2*u^11+3*u^10+3*u^9+3*u^8+u^7+2*u^6+u^5+u^4+u^3+u^2+u+1)*z^3+O(z^4),z,4)

f(p, n) is the coefficient of  of the term u^n*z^p, e.g. f(3, 8) is  the coefficient of u^8*z^3, which is 3.

If you are going to calculate many f(p, n) then you may as well just calculate the series once.

But for a procedure we have

f := proc(p::posint, n::posint)
  local sys, letters, s0, s1, s4, s6, s8, s9, s, w, addition, ser, z, u;
  sys := {letters = 'Union'(s0, s1, s4, s6, s8, s9),
         s0 = 'Atom', s1 = 'Atom', s4 = 'Atom', s6 = 'Atom', s8 = 'Atom', s9 = 'Atom', w = 'Set'(letters)};
  addition := {s(s0) = 0, s(s1) = 1, s(s4) = 4, s(s6) = 6, s(s8) = 8, s(s9) = 9, s(w) = 'Set'(s(letters))};
  Order := p + 1;
  ser := combstruct:-agfseries(sys, addition, 'unlabeled', z, [[u, s]])[w(z, u)];
  ser := collect(convert(ser, 'polynom'), [z,u], 'recursive');
  coeff(coeff(ser, z, p), u, n);
end proc:
  

f(3, 8)

3

CodeTools:-Usage(f(240, 540))

memory used=2.42GiB, alloc change=323.50MiB, cpu time=23.64s, real time=10.31s, gc time=17.11s

1807141

 

 

Download parts.mw

Edit: version that accepts the set of parts as a third argument

f := proc(p::posint, n::{posint, range(posint), list(posint)}, parts::set(nonnegint))
  local sys, letters, s, w, addition, ser, z, u, letts, _s, i, cu;
  letts := seq(cat(_s, i), i in parts);
  sys := {letters = 'Union'(letts), w = 'Set'(letters), letts =~ 'Atom'};
  addition := (map2(apply, s, {letts}) =~ parts) union {s(w) = 'Set'(s(letters))};
  Order := p + 1;
  ser := combstruct:-agfseries(sys, addition, 'unlabeled', z, [[u, s]])[w(z, u)];
  ser := collect(convert(ser, 'polynom'), [z,u], 'recursive');
  cu := coeff(ser, z, p);
  if n::posint then
    coeff(cu, u, n)
  else
    [seq(coeff(cu, u, i), i = n)]
  end if
end proc:

e.g., f(15, 100, {0, 1, 4, 6, 8, 9}); gives 129.

f(15, 90 .. 100, {0, 1, 4, 6, 8, 9}); gives [201, 187, 186, 175, 172, 159, 158, 147, 143, 133, 129]

f(15, [90, 92, 100], {0, 1, 4, 6, 8, 9}); gives [201, 186, 129]

The following is just a programatic version of how you made U. When you say "generalize to higher dimensions" I'm assuming you mean larger matrices, not things in more than two dimensions that wouldn't be matrices.

MMflatten := proc(Q::Matrix(Matrix))
  local m, n, i, j;
  m, n := op(1, Q);
  <seq(`<|>`(seq(Q[i, j], j = 1 .. n)), i = 1 .. m)>
end proc:

 

The first error message says "Warning, expecting only range variable x1 in expression 29239.76608*Ex*x1+794956140.3-29239.76608*Ma-131578.9474*Ey to be plotted but found names [Ex, Ey, Ma]". So you have not given numerical values to Ma, Ex, and Ey. The thing you are trying to plot is

29239.76608*Ex*x1 + 7.949561403*10^8 - 29239.76608*Ma - 131578.9474*Ey

but it can have only the name x1 in. You can try, say 

and you do get a plot, but I don't know what values of these parameters would be sensible.

restart;

eqqq := map(simplify,{sqrt(beta*rho*xi[8]^2 - beta*xi[8]^2)/xi[8] = sqrt(beta*rho - beta), (-rho + 1 + sqrt(beta*rho*xi[8]^2 - beta*xi[8]^2))/(xi[3]*xi[8] - 1) = sqrt(beta*rho - beta), -(-rho*xi[3]*xi[8] + sqrt(beta*rho*xi[8]^2 - beta*xi[8]^2) + xi[3]*xi[8])/(xi[8]*(xi[3]*xi[8] - 1)) = rho - 1}):

In general, simplification on an expression to give zero is easier than comparing both sides of an equation. So I convert the equations to solve to expressions. (I could do this on the solutions after solving.)

eq3:=map(eq->(rhs-lhs)(eq),eqqq);

{(beta*(rho-1))^(1/2)-(beta*xi[8]^2*(rho-1))^(1/2)/xi[8], (beta*(rho-1))^(1/2)-(-rho+1+(beta*xi[8]^2*(rho-1))^(1/2))/(xi[3]*xi[8]-1), rho-1-(-(beta*xi[8]^2*(rho-1))^(1/2)+(rho-1)*xi[3]*xi[8])/(xi[8]*(xi[3]*xi[8]-1))}

sol:=solve(eq3,{xi[8],xi[3]});

{xi[3] = (1+(beta*(rho-1))^(1/2)-rho)/(beta*(rho-1))^(1/2), xi[8] = -(1+(beta*(rho-1))^(1/2)-rho)/(rho-1)}

simplify(subs(sol,eq3));

{(-(beta*(rho-1))^(1/2)+(beta*(1+(beta*(rho-1))^(1/2)-rho)^2/(rho-1))^(1/2)+beta)*(beta*(rho-1))^(1/2)/(rho-(beta*(rho-1))^(1/2)+beta-1), -(-(beta*(rho-1))^(1/2)+(beta*(1+(beta*(rho-1))^(1/2)-rho)^2/(rho-1))^(1/2)+beta)*(rho-1)/(-1-(beta*(rho-1))^(1/2)+rho), -(beta*(rho-1))^(1/2)*(-(beta*(rho-1))^(1/2)+(beta*(1+(beta*(rho-1))^(1/2)-rho)^2/(rho-1))^(1/2)+beta)*(rho-1)/((rho-(beta*(rho-1))^(1/2)+beta-1)*(-1-(beta*(rho-1))^(1/2)+rho))}

We didn't get zero, and it looks like the square roots are a problem. It seems to prevent negatives under the sqrt in sqrt(beta*(rho-1)) we could have

simplify(subs(sol,eq3)) assuming beta>0,rho>1;

{beta^(1/2)*(rho-1)^(1/2)*(1-signum(-1-beta^(1/2)*(rho-1)^(1/2)+rho)), (-1+signum(-1-beta^(1/2)*(rho-1)^(1/2)+rho))*(rho-1)*beta/(beta^(1/2)*(rho-1)^(1/2)-beta-rho+1), (-1+signum(-1-beta^(1/2)*(rho-1)^(1/2)+rho))*(-1-beta^(1/2)*(rho-1)^(1/2)+rho)*beta/(-beta^(1/2)*(rho-1)^(1/2)+beta+rho-1)}

or

simplify(subs(sol,eq3)) assuming beta<0, rho<1;

{0}

If we did want beta>0, rho>1 then the signum is telling us what additional assumption we need

simplify(subs(sol,eq3)) assuming beta>0, rho>1, -1 - sqrt(beta)*sqrt(rho - 1) + rho>0;

{0}

The other way doesn't work

simplify(subs(sol,eq3)) assuming beta>0, rho>1, -1 - sqrt(beta)*sqrt(rho - 1) + rho<0;

{2*beta^(1/2)*(rho-1)^(1/2), 2*beta*(rho-1)/(-beta^(1/2)*(rho-1)^(1/2)+beta+rho-1), -2*beta*(-1-beta^(1/2)*(rho-1)^(1/2)+rho)/(-beta^(1/2)*(rho-1)^(1/2)+beta+rho-1)}

NULL

Download solve.mw

I don't think you can expect sum to give an analytical result, but you can use add to get a numerical result for summing up a large but finite number of terms. (Of course you can use ":" before the evalf to suppress the fraction.) There are probably ways of speeding this up if you want more terms and are only interested in a floating answer.

restart

N := 1000

1000

add(add(ifelse(igcd(m, n) = 1, 1/(m*n)^2, 0), m = 1 .. N), n = 1 .. N); evalf(%)

5532445055376352006682628927753145528970667544410067992788570677846174322152934183442457503057224363774497507108665400532629550147843469142124226144969302173888708851044382140074091091625372693584503837173054948044029528919671403205510054993010871918037249386694096481338979359318665748778806702963161427960903337340452417488905776370472967032563487798981104129939895109772524045706301025651789483034868390590922231566748924949345613510806389477671058222477327056360065465369147452780198772560906625740852672775258980646385235690587185435917247758298832233771101650765956189840592453209825297894908550280663940637781299927315001515221384354297300290305482815046958809041292631478681328638837766548333582204272087678340598645177195671315533413220497463591109093306748119462417424109251101985974126660902042939400696410680608319010223720314403587072423901839/2215401765060199076871263557853118332227324601140256306024979859440747275310142812207631932825650126671317379797411656038492320603818343789239712164205258985121991061929760226402685186979855643381508827056625534669188365623376886098004863359138577635190930503803887889088780257560088147168755709313653758428025570278533530080776685678078392356194935299124340796410605389047851560210858491934810059642988452458045697112532780427976931584727648213025020456341176857124077097214832965005503367731448481531372552150464264617839463298112873781392628562590811184388683738493384943941907659916412512149340044747767907762101305160442358620219612360864537601408693594131025312188123291925844326343537050086070608574762214943201652242835382225022686382300632553803831113114236647288464851072241277028587590660356505813282143777727693670787533236697203835602688000000

2.497264895

NULL

Download igcd.mw

I usually use a colon at the end of an assignment statement storing a plot in a variable, then display it:

plt0 := implicitplot3d(pln, x = -3 .. 3, y = -3 .. 3, z = -3 .. 3, style = patchnogrid, transparency = 0.6):
display(plt0);

What you see with a semicolon in the first line has changed over various versions (plot structure, miniplot, full plot); I didn't check here. But also I would store the plot and not the display of the plot.

That format is hardcoded somewhere hard to find. Here's a hack that takes the existing output and substitutes different strings. Here I use the 5.1f format to get one decimal place.

VennHack.mw

You can collect your Usol values in a table by indexing them. See in red.

ode.mw

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