Advanced Search

dharr

Dr. David Harrington

8145 Reputation

22 Badges

20 years, 326 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

Comments...

dharr 8145
June 18 2025
0 0

@KIRAN SAJJAN Your sahin_base_paper.mw won't run for me - in one of the boundary conditions (the one mentioned in my answer) you have "phis" instead of phi. The plots in that worksheet still have positive H at eta = 0 and you did not change the sign in the boundary condition. lt appears that you didn't read the comment at the end of my worksheet.

Remark...

dharr 8145
June 18 2025
0 0

@mmcdara Thanks. I used to use approxsoln fairly frequently when I was solving these sorts of Navier-Stokes related systems, and it was always hit-and-miss. I think I was lucky here since the actual solution for H at eta=0 was only a small positive value not too different (relatively) from the negative value in the approxsoln. One can provide a digitized version if as in this case you know the solution, but usually that is not the case. Yes, ability to specify just one of the functions would be a nice feature.

Length of output...

dharr 8145
June 18 2025
0 0

@janhardo I agree that this is part of it. Large output should be suppressed because it takes time to decide it is too large to display.

set lambda = 0...

dharr 8145
June 18 2025
0 0

@salim-barzani You omitted the eval(  , lambda=0) part, though when I added that it was still too slow.

plot...

dharr 8145
June 13 2025
0 0

@salim-barzani Here's how to do the plot. I don't know how to export the table; I usually generate data in Maple and write it out to a file for use in some other software. Hope someone else can help you.

 table-e2.mw

exact solution...

dharr 8145
June 12 2025
0 0

@janhardo Thanks for the exact solution and source. I'll work some more on comparisons.

infinite sum...

dharr 8145
June 12 2025
0 0

@salim-barzani You have to recognize the general term manually. For simple cases guessgf can help (see here for an example).

Download Ex1.mw

complex...

dharr 8145
June 12 2025
0 0

@janhardo My procedure should work for cases with conjugates, which is what I think you mean by complex. My only problem is that I don't have any test cases with conjugates that have known exact solutions so I can check it.

b and b[1] etc...

dharr 8145
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 8145
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 8145
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 8145
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 8145
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 8145
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 8145
June 10 2025
0 0

15 minutes sounds about right.

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