dharr

Dr. David Harrington

8330 Reputation

22 Badges

21 years, 3 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 replies submitted by dharr

@astroverted Nice work on the faster method. Use display to combine plots of all kinds.

DQ4.mw

For a 3D plot, I think you will have to use the output=Array option of dsolve, and then use matrixplot or surfdata. I'll try that later when I have a little more time.

The help page for type/float suggests there is support only for double precision. In that case you can do things like in the attached.

Download float.mw

Edit: see my other response for a solution.

@AHSAN Your solution has no variable y. Do you mean x? For phase plots like you show using the commands @tomleslie mentioned, you would need to supply the DE and the initial condions you are interested in.

@mmcdara Sorry if I misunderstood what you meant. I was initially confused about what the OP wanted. My understanding of epsilon is that the contour integration has to be along the real axis but has to include the pole at the origin (the only pole along the real axis in this case). I would write this with integration limits x=-infinity+I*epsilon..+infinity+I*epsilon and leave the epsilon out of the integrand, with the assumption that you need to take the limit as epsilon tends to zero from above. 

Doing it by residues bypasses the details of the contour other than making sure the poles along the real axis are included. I think to do it the hard way would mean having a small semicircle around the poles (into the upper half plane) and taking the limit as their radii went to zero. 

Having said that, I'm still not clear that this is the right definition of the FT, but there seem to be many variations.

Edit: Now I thnk about it, maybe the contour is for the upper half plane excluding the poles along the real axis. That solves the sign problem.

@mmcdara The OP doesn't actually want the Fourier transform of FTI, rather FTI is the definition of the Fourier transform of 1/sinh(x)^2.

@mmcdara I liked your hand solution and voted it up. But since @yangtheary never seems to respond, I agree that it may not be worthwhile.

@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:
DQ.mw

@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.
 

NULL

restart

with(DifferentialGeometry):
NULL

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

@WA573 
 

restart

with(PDEtools):

declare(u(x, t));

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

(1)

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

(2)

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

(3)

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

(4)

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

(5)

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

Note that

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

true

(6)

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

(7)

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

(8)

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.

First 53 54 55 56 57 58 59 Last Page 55 of 87