Dr. David Harrington

3481 Reputation

17 Badges

18 years, 93 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 professor of chemistry at the University of Victoria, BC, Canada, where 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

@astroverted Once I fixed the syntax errors, Maple told me things like "expected 1 boundary condition but found 2" until I had only one initial condition and one boundary condition. That's what I expected since there was one equation first order in time and one first order in z. The other "initial" and "boundary" conditions in the paper are simply found from the others so are not independent. Eq (16) is just putting M=1 (at t=0) into Eq (13) and solving the ode; Eq (18) is from putting I=I0 into Eq (14) and solving.

Maple can't solve equations of this form as it only solves time-based systems. But it can be solved, for example by the method of lines or differential quadrature. Here is an example:

@zenterix How is that person executing the .mpl file? By running command line Maple? Or (as you suggest above) by reading it into a worksheet? The way I generate packages or help databases is by running a worksheet which gets other files as necessary from the same directory, and uses worksheetdir. (Getting them from a subdirectory is a simple modification of that.) Someone else can run that worksheet which is in the same directory as all the required files on their system, since I would give then the zipped contents of the directory.

If you read a .mpl file with read, then the code is running in the worksheet you read it into, not in the .mpl file.

@mmcdara Yes, I was developing some code like that before I realized I was reinventing the wheel. You can use generator=(()->1) so you always get 1, then density gives the probability. Actually not quite because the ones along the diagonal were removed. Since the native graph datastructure has to be made from the matrix, this may not be as efficient as directly generating the edges and vertices.

@zenterix That's why I'm confused. You are opening a worksheet to read the .mpl file. So in order to read the .mpl file at this point you know where it is. So what do you mean by getting the path "programatically from within the mpl file"?

Note that once the worksheet is saved and then you exit Maple and reopen, the worksheet is at currentdir(). So if I distribute some code in a series of files, I can distribute a worksheet in the same directory that reads test.mpl or does something else. Then when the user runs that worksheet, currentdir() (or worksheetdir) will be where that worksheet is and where test.mpl is.

@JohannesC My Maple version (2015) had some errors so I cut and pasted CONV_2 to a new worksheet. But now you have inserted the missing line, I am getting the same result as you. 

The problem is that rho has an internal form that is different from the way it looks (try lprint(rho) to see this). The workaround is to write in the expression for rho in the line that generates CONV_2 - then simplify works.




DGsetup([r, theta, phi], M)

`frame name: M`

Defining the Metric Tensor and Calculating the Connection Coefficients


g := evalDG(`&t`(dr, dr)+r^2*(sin(phi)^2)*`&t`(dtheta, dtheta)+r^2*`&t`(dphi, dphi))

_DG([["tensor", M, [["cov_bas", "cov_bas"], []]], [[[1, 1], 1], [[2, 2], r^2*sin(phi)^2], [[3, 3], r^2]]])

g_inv := InverseMetric(g)

_DG([["tensor", M, [["con_bas", "con_bas"], []]], [[[1, 1], 1], [[2, 2], 1/(r^2*sin(phi)^2)], [[3, 3], 1/r^2]]])

C := Christoffel(g):

DifferentialGeometry:-Tools:-DGinfo(C, "ObjectComponents")

[[[1, 2, 2], -r*sin(phi)^2], [[1, 3, 3], -r], [[2, 1, 2], 1/r], [[2, 2, 1], 1/r], [[2, 2, 3], cos(phi)/sin(phi)], [[2, 3, 2], cos(phi)/sin(phi)], [[3, 1, 3], 1/r], [[3, 2, 2], -sin(phi)*cos(phi)], [[3, 3, 1], 1/r]]

Defining the PDF as a Tensor (Metric) Density


f := f__1(r, theta, phi)

f__1(r, theta, phi)

rho := MetricDensity(g, -1)

_DG([["tensor", M, [[], [["bas", -1]]]], [[[], 1/(r^4*sin(phi)^2)^(1/2)]]])

Tools:-DGinfo(rho, "TensorDensityType")

[["bas", -1]]

Convective Operator


v := DifferentialGeometry:-evalDG(VectorCalculus:-`+`(VectorCalculus:-`+`(VectorCalculus:-`*`(v__r, D_r), VectorCalculus:-`*`(v__theta, D_theta)), VectorCalculus:-`*`(v__phi, D_phi)))

_DG([["vector", M, []], [[[1], v__r], [[2], v__theta], [[3], v__phi]]])

CONV := CovariantDerivative(VectorCalculus:-`*`(v, f), C):

CONV_1 := simplify(ContractIndices(CONV, [[1, 2]])):

CONV_2 := expand(eval(CONV_1, f__1(r, theta, phi) = VectorCalculus:-`*`(f__2(r, theta, phi), 1/sqrt(VectorCalculus:-`*`(r^4, sin(phi)^2)))))

cos(phi)*f__2(r, theta, phi)*v__phi/(sin(phi)*(r^4*sin(phi)^2)^(1/2))+v__phi*(diff(f__2(r, theta, phi), phi))/(r^4*sin(phi)^2)^(1/2)-sin(phi)*r^4*v__phi*f__2(r, theta, phi)*cos(phi)/(r^4*sin(phi)^2)^(3/2)+(diff(f__2(r, theta, phi), theta))*v__theta/(r^4*sin(phi)^2)^(1/2)+v__r*(diff(f__2(r, theta, phi), r))/(r^4*sin(phi)^2)^(1/2)-2*sin(phi)^2*r^3*v__r*f__2(r, theta, phi)/(r^4*sin(phi)^2)^(3/2)+2*f__2(r, theta, phi)*v__r/(r*(r^4*sin(phi)^2)^(1/2))

CONV_2 := simplify(CONV_2)

-(cos(phi)^2-1)*(v__phi*(diff(f__2(r, theta, phi), phi))+(diff(f__2(r, theta, phi), theta))*v__theta+v__r*(diff(f__2(r, theta, phi), r)))/(sin(phi)^2*(r^4*sin(phi)^2)^(1/2))



Download Tensor_Calc_simplified2.mw




declare(u(x, t));

u(x, t)*`will now be displayed as`*u


PDE := diff(u(x, t), t)-(diff(u(x, t), `$`(x, 2), t))+3*u(x, t)^2*(diff(u(x, t), x))-2*(diff(u(x, t), x))*(diff(u(x, t), `$`(x, 2)))-u(x, t)*(diff(u(x, t), `$`(x, 3)));

diff(u(x, t), t)-(diff(diff(diff(u(x, t), t), x), x))+3*u(x, t)^2*(diff(u(x, t), x))-2*(diff(u(x, t), x))*(diff(diff(u(x, t), x), x))-u(x, t)*(diff(diff(diff(u(x, t), x), x), x))


ODE := convert(simplify(eval(PDE, u(x, t) = U(t*lambda__2+x*lambda__1)), {t*lambda__2+x*lambda__1 = eta}), diff);

(-2*(diff(U(eta), eta))*(diff(diff(U(eta), eta), eta))-(diff(diff(diff(U(eta), eta), eta), eta))*U(eta))*lambda__1^3-(diff(diff(diff(U(eta), eta), eta), eta))*lambda__1^2*lambda__2+3*U(eta)^2*(diff(U(eta), eta))*lambda__1+(diff(U(eta), eta))*lambda__2


ODE1 := simplify((-2*(diff(U(eta), eta))*(diff(diff(U(eta), eta), eta))-(diff(diff(diff(U(eta), eta), eta), eta))*U(eta))*lambda__1^3-(diff(diff(diff(U(eta), eta), eta), eta))*lambda__1^2*lambda__2+3*U(eta)^2*(diff(U(eta), eta))*lambda__1+(diff(U(eta), eta))*lambda__2, 'symbolic')

-2*(diff(U(eta), eta))*(diff(diff(U(eta), eta), eta))*lambda__1^3-(diff(diff(diff(U(eta), eta), eta), eta))*U(eta)*lambda__1^3+3*U(eta)^2*(diff(U(eta), eta))*lambda__1-(diff(diff(diff(U(eta), eta), eta), eta))*lambda__1^2*lambda__2+(diff(U(eta), eta))*lambda__2


anst := U(eta) = sum(a[i]*csc(eta)^i, i = -2 .. 2);

U(eta) = a[-2]/csc(eta)^2+a[-1]/csc(eta)+a[0]+a[1]*csc(eta)+a[2]*csc(eta)^2


Eq := normal(subs(anst, ODE1)):

Note that

is(cot(x)^2 = csc(x)^2-1);



So we can convert the cot to csc and leave only powers of csc(eta) and lambda with

Eq2 := simplify(numer(Eq)/cot(eta), {cot(eta)^2 = csc(eta)^2-1}):

vars := `minus`(indets(Eq2), {eta, csc(eta)});

{lambda__1, lambda__2, a[-2], a[-1], a[0], a[1], a[2]}


Write csc(eta) as x

Eq3 := eval(Eq2, csc(eta) = x):

solve(identity(Eq3, x), vars);

{lambda__1 = lambda__1, lambda__2 = lambda__2, a[-2] = 0, a[-1] = 0, a[0] = a[0], a[1] = 0, a[2] = 0}, {lambda__1 = lambda__1, lambda__2 = (-3/2+(1/2)*(-128*lambda__1^4+9)^(1/2))*lambda__1, a[-2] = 0, a[-1] = 0, a[0] = -(8/3)*lambda__1^2-1/2+(1/6)*(-128*lambda__1^4+9)^(1/2), a[1] = 0, a[2] = 8*lambda__1^2}, {lambda__1 = lambda__1, lambda__2 = (-3/2-(1/2)*(-128*lambda__1^4+9)^(1/2))*lambda__1, a[-2] = 0, a[-1] = 0, a[0] = -(8/3)*lambda__1^2-1/2-(1/6)*(-128*lambda__1^4+9)^(1/2), a[1] = 0, a[2] = 8*lambda__1^2}, {lambda__1 = 0, lambda__2 = 0, a[-2] = a[-2], a[-1] = a[-1], a[0] = a[0], a[1] = a[1], a[2] = a[2]}


Download CM_TW2.mw

@zenterix I added an example to my answer, using LibraryTools as @acer pointed out.

@AHSAN So you will need to arrange your loops so that the value of Q and lambda are known at the time you call dsolve.

An outer loop can change values of Q and lambda and the innner loop with the other parameters as you have now.

@AHSAN Suggest you fix it then.

You define the function theta(y) in terms of y. Then you use dsolve to find u(y), but you give boundary conditions for the known function theta(y), not the unkown function u(y).

@moh111 What do you know about the signs or ranges of k, re and the variables?

@tomleslie The OP specified the 8 unknowns they were interested in. My understanding is that the OP wants these as functions of k and re. Setting infolevel[solve]:=5 shows that solve uses methods for polynomial systems, but they are very slow. (If k and re are specified then fsolve easily gives a solution.)

You should not use both N(x,t) and N[0] together - use N__0 instead of N[0] (looks the same), and the same for omega, E etc.

You set a Matrix equal to scalar zero; I assume that was meant to be the determinant of that Matrix.

If I guess at what substitution you want, then I get the attached, which still has cosines and x and t in. You will need to specify exactly the steps that lead to the result you expect. My guess is that it is not a matter of substitution, but of solving a PDE and getting some eigenvalues for E, but the equations aren't familiar to me, so many more details are required.


@bstuan I was just too lazy to write out the 2023 case :-)

Edit: That is, (A.B)^n is not equal to A^n.B^n in general, though there are some cases where it is true e.g., commuting matrices, n=1, n=0.

@bstuan Note (A.B)^2 = A.B.A.B <> A.A.B.B

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