Advanced Search

dharr

Dr. David Harrington

8205 Reputation

22 Badges

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

Social Networks and Content at Maplesoft.com

Maple Application Center
Contact Dr. David Harrington
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 replies submitted by dharr

b and b[1] etc...

dharr 8205
June 11 2025
0 0

@salim-barzani You can't have both indexed b (b[1], b[2], b[3]) and unindexed b together as parameters. Change b to b[0]. If you can get this allegedly exact solution to satisfy pdetest (for any parameter set and/or x and t values) then I will look more closely at it.

only U(x,t) and conjugate(U(x,t))...

dharr 8205
June 11 2025
0 0

@janhardo Your modification expands only in U(x,t) and conjugate(U(x,t)) = V(x,t), but the papers also use the derivatives of these wrt x. Even without conjugates, the example with the code expands in both U(x,t) and diff(U(x,t),x), and gives a correct solution.

problem not set up correctly...

dharr 8205
June 11 2025
0 0

 I tried the code, but the initial condition doesn't agree with Fig. 1 so something is wrong before the code runs.

Laplace_Adomiian_Decomposition_procedure4.mw

preliminary code...

dharr 8205
June 10 2025
0 0

@dharr Here is a preliminary procedure that finds an exact solution, but for a simple case without abs or conjugates. The code to find vars is fragile, since I'm currently not sure what cases can be done by this method, so use infolevel[LAD] := 1 or 2 and look at the expansion variables to check it makes sense.

Laplace Adomian Decomposition method. Routine LAD is in the startup region - code reproduced below.

> 

restart

Recommended to load Physics package if derivatives contain abs or conjugate.

> 

NULL

PDE in diff notation - Sec 4, Wazwaz paper

> 

pde := diff(U(x, t), t)+U(x, t)^2*(diff(U(x, t), x))

diff(U(x, t), t)+U(x, t)^2*(diff(U(x, t), x))

> 

infolevel[LAD] := 2

LAD takes a pde, unknown function, initial condition, time variable, and degree of approximation (N=10 means the approximate solution is add(U[i](x,t), i=0..N)

> 

soln := LAD(pde, U(x, t), 3*x, t, 10)

LAD: L = diff(U(x,t),t), R = 0, N = -U(x,t)^2*diff(U(x,t),x)

LAD: expansion variables are: [U(x,t), diff(U(x,t),x)]

175692092397588*t^10*x^11-5650915252554*t^9*x^10+184670433090*t^8*x^9-6155681103*t^7*x^8+210450636*t^6*x^7-7440174*t^5*x^6+275562*t^4*x^5-10935*t^3*x^4+486*t^2*x^3-27*t*x^2+3*x

Extrapolate to an infinite number of terms.

> 

algsubs(x*t = y, expand(soln*t)); ListTools:-Reverse([coeffs(%, y), 0]); approx := subs(y = x*t, gfun:-guessgf(%, y)[1])/t

175692092397588*y^11-5650915252554*y^10+184670433090*y^9-6155681103*y^8+210450636*y^7-7440174*y^6+275562*y^5-10935*y^4+486*y^3-27*y^2+3*y

[0, 3, -27, 486, -10935, 275562, -7440174, 210450636, -6155681103, 184670433090, -5650915252554, 175692092397588]

6*x/(1+(36*t*x+1)^(1/2))

This actually is an exact solution

> 

simplify(eval(pde, U(x, t) = approx))

0

> 

LAD := proc(pde, Uxt::function(name), inx, t::name, N::posint)
option `D.A. Harrington Jun 2025 v 1.03`;
local Ru, u, pde1, Ut, q1, q2, scale, term, T, RRu, NNu, A, s, lambda, vars, nvars, vvars, NNv, i, v, vals, xt, n, j, Ai, soln, z;
# "with(Physics)" recommended before calling if derivatives contain abs or conjugate
if has(inx, t) then error "initial condition cannot contain %1", t end if;
if type(pde1, `=`) then pde1 := (lhs-rhs)(expand(pde)) else pde1 := expand(pde) end if;
if not type(pde, `+`) then error "pde expected to be a sum of terms" end if;
Ut := diff(Uxt,t);        # currently higher derivatives wrt t not handled
q1, q2 := selectremove(has, pde1, Ut); # q1 can be any derivative wrt t
if q1 = 0 then error "derivative of %1 wrt %2 not found", Uxt, t end if;
if indets((scale := q1/Ut)) <> {} then error "time-dependent term can only be %1 multiplied by a constant", Ut end if;
if (q1 := remove(has, q2 + z, Uxt) - z) <> 0 then error "homogeneous terms %1 not handled", q1 end if;
q2 := expand(-q2/scale);
RRu, NNu := selectremove(term -> not has(eval(term, Uxt = T*Uxt)/T, T), q2 + z);        # linear and nonlinear parts; add z to force type `+`
NNu := NNu - z;                # correct for +z
if NNu = 0 then error "pde has no nonlinear part" end if;
userinfo(2, procname, "L = %1, R = %2, N = %3", Ut, RRu, NNu);
xt := op(Uxt);
u[0] := inx;
soln := inx;
vars := select(has, [op(indets(NNu))], Uxt);
userinfo(1, procname, "expansion variables are: %1", vars);
nvars := nops(vars);
vvars := [seq(v[i], i = 1 .. nvars)];
NNv := subs(vars =~ vvars, NNu);
vals[0] := simplify(eval(subs(Uxt = u[0], vars))) assuming real;
Ai := eval(subs(vvars =~ vals[0], NNv));
A[0] := simplify(Ai) assuming real;
Ru := simplify(eval(subs(Uxt = u[0], RRu))) assuming real;
for n to N do
        u[n] := inttrans:-invlaplace(inttrans:-laplace(Ru + A[n - 1], t, s)/s, s, t);
        soln := soln + u[n];
        vals[n] := simplify(eval(subs(Uxt = u[n], vars))) assuming 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)/n!, lambda = 0));
        Ru := simplify(eval(subs(Uxt = u[n], RRu))) assuming real;
end do;
soln
end proc:
> 

NULL

Download Laplace_Adomiian_Decomposition_procedure3b.mw

not exact?...

dharr 8205
June 10 2025
0 0

@salim-barzani Eqn 9 does not seem to be an exact solution - significant error for small x and t. If you have one that is exact for which the LADM method can be tested against, I will work on a general procedure for the LADM method.

> 

restart
> 

with(inttrans)
> 

with(PDEtools)
> 

with(DEtools)
> 

with(Physics)
> 

declare(u(x, t), quiet); declare(v(x, t), quiet); declare(U(x, t), quiet)
> 

undeclare(prime)

`There is no more prime differentiation variable; all derivatives will be displayed as indexed functions`

> 

pde0 := I*(diff(U(x, t), t))+I*a[1]*(diff(U(x, t), `$`(x, 3)))+a[2]*(diff(U(x, t), `$`(x, 4)))+U(x, t)*conjugate(U(x, t))*(b*U(x, t)+I*sigma*(diff(U(x, t), x)))-I*(alpha*(diff(U(x, t), x))+beta*(diff(U(x, t)^2*conjugate(U(x, t)), x))+mu*(diff(U(x, t)*conjugate(U(x, t)), x))*U(x, t))

I*(diff(U(x, t), t))+I*a[1]*(diff(diff(diff(U(x, t), x), x), x))+a[2]*(diff(diff(diff(diff(U(x, t), x), x), x), x))+U(x, t)*conjugate(U(x, t))*(b*U(x, t)+I*sigma*(diff(U(x, t), x)))-I*(alpha*(diff(U(x, t), x))+beta*(2*U(x, t)*conjugate(U(x, t))*(diff(U(x, t), x))+U(x, t)^2*(diff(conjugate(U(x, t)), x)))+mu*((diff(U(x, t), x))*conjugate(U(x, t))+U(x, t)*(diff(conjugate(U(x, t)), x)))*U(x, t))

> 

NULL
> 

A := (-6*k^2*a[2]+3*k*a[1])*sqrt(-30*a[2]/(-beta*k+k*sigma+b))/(10*a[2]); v := 4*k^3*a[2]-3*k^2*a[1]-alpha; B := (1/2)*sqrt(-(-6*k^2*a[2]+3*k*a[1])/(5*a[2])); sol_exact := A*sech(B*(-t*v+x))^2*exp(I*(-k*x+t*w+theta))

(1/10)*(-6*k^2*a[2]+3*k*a[1])*(-30*a[2]/(-beta*k+k*sigma+b))^(1/2)/a[2]

4*k^3*a[2]-3*k^2*a[1]-alpha

(1/10)*(-5*(-6*k^2*a[2]+3*k*a[1])/a[2])^(1/2)

(1/10)*(-6*k^2*a[2]+3*k*a[1])*(-30*a[2]/(-beta*k+k*sigma+b))^(1/2)*sech((1/10)*(-5*(-6*k^2*a[2]+3*k*a[1])/a[2])^(1/2)*(x-(4*k^3*a[2]-3*k^2*a[1]-alpha)*t))^2*exp(I*(-k*x+t*w+theta))/a[2]

> 

params := convert({alpha = 0.8e-1, b = 4.9, beta = -1.4, k = 1, mu = 1.02, sigma = 0.1e-1, theta = 0, w = 0, a[1] = 2.20, a[2] = -.35}, rational)

{alpha = 2/25, b = 49/10, beta = -7/5, k = 1, mu = 51/50, sigma = 1/100, theta = 0, w = 0, a[1] = 11/5, a[2] = -7/20}

> 

q1 := `assuming`([simplify(eval(eval(pde0, params), U(x, t) = eval(sol_exact, params)))], [real]); evalf[50](eval(q1, {t = 2*(1/100), x = 1/1000})); evalf[10](eval(q1, {t = 2, x = 1})); evalf(eval(q1, {t = 17, x = 13}))

(272484/27054125)*exp(-I*x)*(-783293/12528+I*(cosh((1/1750)*(25*x+202*t)*6090^(1/2))^4-3*cosh((1/1750)*(25*x+202*t)*6090^(1/2))^2+6293/5048)*42^(1/2)*tanh((1/1750)*(25*x+202*t)*6090^(1/2))*145^(1/2)*sech((1/1750)*(25*x+202*t)*6090^(1/2))^4)*631^(1/2)*42^(1/2)*sech((1/1750)*(25*x+202*t)*6090^(1/2))^2

-99.235569661618325350229067688990397347697093307085-16.239718062808372129333287933226454261993028399355*I

0.5055914512e-14+0.1502809341e-13*I

-0.3957818373e-143+0.1604549128e-142*I

> 

NULL

Download test_exact_soln.mw

pdetest...

dharr 8205
June 10 2025
0 0

@salim-barzani I think you need a paper with an exact solution that you can verify by pdetest that it is correct. Is that the case for the solution here?

yes...

dharr 8205
June 10 2025
0 0

15 minutes sounds about right.

more progress?...

dharr 8205
June 09 2025
0 0

@dharr The paper keeps 4 variables in A[i], with both diff(U(x,t),x) and diff(conjugate(U(x,t)),x) present. When I do that, I get the wrong coefficients for A[2], but the right number of terms, and the sums of subcripts = 2 as required. On the other hand in your post here you get the coefficients correct for A[2], but you do not have any diff(conjugate(U(x,t)),x).

There is definitely some difficulty about processing diff(conjugate(U(x,t)),x) that I do not understand - I have been careful to substitute for U(x,t) then conjugate so that the differentiation sees only expressions without conjugate. It's not clear to me how to proceed from here.

The whole process with the lambda's can be replaced by some recurrence formulas, but these have partial differentiation wrt the v[i] and again the conjugate variables make this hard to work out.

Laplace Adomian Decomposition method.

> 

restart
> 

with(Physics); with(inttrans)
> 

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]); 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)]

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 seems to agree with A[1] below Eq 25 but I didn't check every term.
However 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(v[j][i]*lambda^i, i = 0 .. n), j = 1 .. nvars)}, NNv), `$`(lambda, 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)), 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)*v[1][1]*v[2][0]*v[4][0]+(2*I)*v[1][0]*v[2][1]*v[4][0]+(2*I)*v[1][0]*v[2][0]*v[4][1]+(4*I)*v[1][0]*v[3][0]*v[1][1]+(2*I)*v[1][0]^2*v[3][1]+(3*I)*v[1][0]^2*v[4][0]^2*v[1][1]+(2*I)*v[1][0]^3*v[4][0]*v[4][1]

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

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

Next iteration; of course this can be put in a for loop.

> 

n := 2; 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(v[j][i]*lambda^i, i = 0 .. n), j = 1 .. nvars)}, NNv), `$`(lambda, 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)), lambda = 0)); Ru[n] := `assuming`([simplify(eval(RRu, U(x, t) = u[n]))], [real])

2

-(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]

(12*I)*v[1][0]^2*v[4][0]*v[1][1]*v[4][1]+(4*I)*v[1][0]*v[2][0]*v[4][2]+(4*I)*v[1][1]*v[2][0]*v[4][1]+(6*I)*v[1][0]^2*v[4][0]^2*v[1][2]+(4*I)*v[1][0]^2*v[3][2]+(8*I)*v[1][0]*v[3][0]*v[1][2]+(8*I)*v[1][0]*v[3][1]*v[1][1]+(4*I)*v[1][0]*v[2][2]*v[4][0]+(4*I)*v[1][0]^3*v[4][0]*v[4][2]+(4*I)*v[1][2]*v[2][0]*v[4][0]+(2*I)*v[1][0]^3*v[4][1]^2+(4*I)*v[1][0]*v[2][1]*v[4][1]+(4*I)*v[1][1]^2*v[3][0]+(4*I)*v[1][1]*v[2][1]*v[4][0]+(6*I)*v[1][0]*v[4][0]^2*v[1][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

 

 

Download Laplace_Adomiian_Decomposition_method_3.mw

Progress...

dharr 8205
June 06 2025
0 0

@dharr This is how I would approach it, but there is something I'm not understanding about Eq 8.12 (aside from it suggesting setting lambda=0 before differentiating wrt lambda. I'm out of time for the next few days.

Laplace Adomian Decomposition method.

> 

restart
> 

with(Physics); with(inttrans)
> 

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(U(x, t), x, x)+2*(Diff(U(x, t)*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(U(x, t), x, x))

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-I*(Diff(U(x, t), x, x))

u[0] starts as the initial condition.

> 

u[0] := beta*exp(I*x); A[0] := `assuming`([simplify(value(eval(NNu, U(x, t) = u[0])))], [real])

beta*exp(I*x)

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

u[1] is a function of t for the next Laplace transform.

> 

Ru[0] := `assuming`([simplify(value(eval(RRu, U(x, t) = u[0])))], [real]); u[1] := invlaplace(laplace(Ru[0]+A[0], t, s)/s, s, t)

-I*beta*exp(I*x)

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

Using Eq 8.12 - something wrong here or with the formula?

> 

n := 1; `assuming`([(eval(diff(value(eval(NNu, U(x, t) = add(lambda^i*u[i], i = 0 .. n))), `$`(lambda, n)), lambda = 0))/factorial(n)], [real]); A[n] := `assuming`([simplify(eval(%))], [real])

1

(3*I)*u[0]^2*conjugate(u[0])^2*u[1]+(2*I)*u[0]^3*conjugate(u[0])*conjugate(u[1])

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

> 

NULL
> 

NULL

Download Laplace_Adomiian_Decomposition_method_2.mw

uploading...

dharr 8205
June 06 2025
0 0

@rlewis Click the green up-arrow; then choose the file, click upload, then click insert link. (Sometimes insert content doesn't work, though the link will usually work.

error...

dharr 8205
June 05 2025
0 0

@janhardo You have 

uxx := diff(u[0](x), x $ 2);

but u[0] is already a function of x, so it should be

uxx := diff(u[0], x $ 2);

This is the basic way to do it, with everything in a big loop.

function of t...

dharr 8205
June 05 2025
0 0

@salim-barzani As I said before, 

-diff(u[i], x $ 2)*I + A[i]

should have some t dependence or you can't expect to correctly take a Laplce transform wrt t. So you need to put the dependence of t in here (or just have it in from the beginning).

P1...

dharr 8205
June 05 2025
0 0

@salim-barzani I didn't see where you divided by 2. P[1] looks the same in both versions?

Download P1.mw

Workaround...

dharr 8205
June 05 2025
0 0 
int(convert(1/A, Heaviside), x=-1..1);

Bad T[n]...

dharr 8205
June 05 2025
0 0

@salim-barzani Your T[n] in DrD.mw has conjugate(u[k]) outside the summation, so it just stays as k.

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