Dr. David Harrington

2403 Reputation

15 Badges

17 years, 146 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 answers submitted by dharr

Well, you are a moderator (note the symbol below your avatar), "a selected group of long-time MaplePrimes users with a reputation for positive contributions." - see https://www.mapleprimes.com/help/moderation.

There are some good reasons to edit others' posts; removing copyrighted content is one I've used. Another is removing duplicate copies of posted worksheets.

There are guidelines on the cited page. One of them is:

"In general, try to avoid making spelling or grammar changes unless they are particularly egregious, but feel free to edit message formatting if it messes up the display of the page." Since I don't think @acer's grammar was egregious, perhaps you are in violation :)

But seriously, this policy does not seem to have led to much difficulty.

@dharr OK, I see now what you want - here is how I would do it with laplace. Not sure about with Elziki



solve(invEqs, {theta[1], theta[2], theta[3]}, explicit) gives two solutions, though they are many pages long. Once you have theta[1], theta[2] and theta[3], you can find sin(theta[1]) or various other derived quantities.

Your 7x7 Matrix has only Rank 6, so it must have Determinant zero, no matter what the values of k, Bi, omega etc are. If this isn't what your expect, you will need to check the construction of the Matrix. I didn't see any obvious reason for this, like two rows or columns the same.

Maple does have this as a builtin command.

CurveFitting:-PolynomialInterpolation(lxi, Lu, xi)

returns -(1/3)*xi^2+(4/3)*xi+1


The dsolve command solves DAE systems. See the help ?dsolve,dae. There are five DEA methods,  rkf45_dae, ck45_dae, rosenbrock_dae, or mebdfi, and if you specify stiff=true, you will get the rosenbrock method. There is a brief desciption and some options described on that help page and the ?dsolve,dae_extension help page.



y1 := exp(x)


y2 := r*x


When the two curves are tangent, then there will be only one solution, which occurs for r=e

solve({y1 = y2, diff(y1, x) = diff(y2, x)}, {r, x});

{r = exp(1), x = 1}

r=e, just touching

plot(eval([y1, y2], r = exp(1)), x = 0 .. 3);

r>e, two solutions

plot(eval([y1, y2], r = exp(1)+.5), x = 0 .. 3);

r<e, no real solutions

plot(eval([y1, y2], r = exp(1)-.5), x = 0 .. 3);




Download Number_solutions.mw

Here is one way. You may have to optimize plot parameters for the other plots.

Download sim_plot.mw

pdsolve can give a general solution, which can probably be worked up with BCs and IC in some cases to give an analytical solution, though not immediately in my version for the BCs you gave (incompatible?) :


interface(version); with(VectorCalculus)

`Standard Worksheet Interface, Maple 2017.3, Windows 8.1, September 27 2017 Build ID 1265877`

pde := diff(u(t, x, y, z), t) = k*Laplacian(u(t, x, y, z), [x, y, z])+l*u(t, x, y, z)+m

diff(u(t, x, y, z), t) = k*(diff(diff(u(t, x, y, z), x), x)+diff(diff(u(t, x, y, z), y), y)+diff(diff(u(t, x, y, z), z), z))+l*u(t, x, y, z)+m


PDESolStruc(u(t, x, y, z) = _F1(t)*_F2(x)*_F3(y)*_F4(z)-m*(_C1+_C2*sin(l^(1/2)*x/k^(1/2))+_C3*cos(l^(1/2)*x/k^(1/2)))/(l*_C1), [{diff(_F1(t), t) = _c[1]*_F1(t), diff(diff(_F2(x), x), x) = _c[2]*_F2(x), diff(diff(_F3(y), y), y) = _c[3]*_F3(y), diff(diff(_F4(z), z), z) = _F4(z)*_c[1]/k-_F4(z)*_c[2]-_F4(z)*_c[3]-_F4(z)*l/k}])

bcs := u(t, 0, y, z) = 0, u(t, L, y, z) = 0, u(t, x, 0, z) = 0, u(t, x, L, z) = 0, u(t, x, y, 0) = 0, u(t, x, y, L) = 0

u(t, 0, y, z) = 0, u(t, L, y, z) = 0, u(t, x, 0, z) = 0, u(t, x, L, z) = 0, u(t, x, y, 0) = 0, u(t, x, y, L) = 0

pdsolve({bcs, pde})



Download pdsolve.mw

Numerical solutions are also possible, see ?pdsolve,numeric. [Edit - only for time and two space variables.]

To make b^1 a name, select it with the mouse, and then use the right-click (ctrl-click on mac)  menu to choose 2-D math -> convert to -> atomic variable.



You can set up a Vector of r values, say rvals, at which you want the solution, and then use output=rvals. The output is a 2x1 Matrix, from which you can extract the information you want.


The option axes=none removes the axes and labels.

If you substitute r=3/2 in Maple's general solution you get zero, so it will never match the initial condition. So Maple's general solution is not general enough; a bug IMO. Going through the standard separation of variables steps leads to an analytical solution.

restart; with(PDEtools); heat := diff(u(r, t), t) = (diff(r^2*(diff(u(r, t), r)), r))/r^2; bc1 := u(1, t) = 0; bc2 := u(2, t) = 0; bc := bc1, bc2; ic := u(r, 0) = -sin(Pi*r); sol := pdsolve([heat, bc], HINT = `*`)

diff(u(r, t), t) = (2*r*(diff(u(r, t), r))+r^2*(diff(diff(u(r, t), r), r)))/r^2

u(1, t) = 0

u(2, t) = 0

u(1, t) = 0, u(2, t) = 0

u(r, 0) = -sin(Pi*r)

u(r, t) = Sum(I*_C1(n)*sin(2*Pi*n*r)*exp(-4*Pi^2*n^2*t)/r, n = 1 .. infinity)

The "general" solution above is always zero at r=3/2

simplify(eval(sol, r = 3/2))

u(3/2, t) = 0

Look for a separation of variable solution

heat2 := simplify((eval(heat, u(r, t) = f(t)*g(r)))/(f(t)*g(r)))

(diff(f(t), t))/f(t) = ((diff(diff(g(r), r), r))*r+2*(diff(g(r), r)))/(r*g(r))

lhs function of t only and rhs function of r only, so each are equal to the separation constant, which we write as -lambda^2

der := -lambda^2 = rhs(heat2)

-lambda^2 = ((diff(diff(g(r), r), r))*r+2*(diff(g(r), r)))/(r*g(r))

If given both boundary conditions, dsolve gives the trivial solution and two others; the third one works up to Maple's general solution above.

sol1 := dsolve({der, g(1) = 0, g(2) = 0}, {lambda, g(r)})

{lambda = _C1, g(r) = 0}, {lambda = 2*Pi*_Z3+Pi, g(r) = _C2*sin((2*Pi*_Z3+Pi)*r)/r}, {lambda = 2*Pi*_Z2, g(r) = _C2*sin(2*Pi*_Z2*r)/r}

Taking both solutions, we have lambda a positive integer, with

gr := a__n*sin(n*Pi*r)/r; leq := lambda = n*Pi


lambda = Pi*n

Use this lambda in the time part

eval(-lambda^2 = lhs(heat2), leq); ft := rhs(dsolve({%, f(0) = 1}))

-Pi^2*n^2 = (diff(f(t), t))/f(t)


And the general solution is the sum over positive n of:

sol4 := ft*gr


So we now want to solve the initial condition part.
The a__n values are different for each term

sol5 := rhs(ic) = eval(Sum(sol4, n = 1 .. infinity), t = 0)

-sin(Pi*r) = Sum(a__n*sin(Pi*n*r)/r, n = 1 .. infinity)

Since we know the sin functions are orthogonal for different n we can see that multiplying by "r*sin((Pi*m*r)" and integrating will send all terms except m=n to zero

`assuming`([int(r*sin(m*Pi*r)*(sin(n*Pi*r)/r), r = 1 .. 2)], [n::posint, m::posint, m <> n])


And for m=n

`assuming`([int(r*sin(n*Pi*r)*(sin(n*Pi*r)/r), r = 1 .. 2)], [n::posint])


 We need to solve

eq := `assuming`([int(r*sin(n*Pi*r)*rhs(ic), r = 1 .. 2) = (1/2)*a__n], [n::posint])

-2*n*((-1)^n+1)/(Pi^2*(n^4-2*n^2+1)) = (1/2)*a__n

an := simplify(solve(eq, a__n))


We see that odd n gives zero and n=1 needs to be treated separately

a1 := solve(int(r*sin(Pi*r)*rhs(ic), r = 1 .. 2) = (1/2)*a__1, a__1)


Inserting nto the general solution gives

sol6 := eval(sol4, {a__n = a1, n = 1})+Sum(eval(sol4, a__n = an), n = 2 .. infinity)

-(3/2)*exp(-Pi^2*t)*sin(Pi*r)/r+Sum(-4*exp(-Pi^2*n^2*t)*n*((-1)^n+1)*sin(Pi*n*r)/(Pi^2*(n-1)^2*(n+1)^2*r), n = 2 .. infinity)

Check it with pdetest - won't simplify to zero so try select r and t values - last value not zero maybe still something not quite right or a numerical error?

evalf(eval(pdetest(u(r, t) = sol6, [heat, bc, ic]), {r = 1.2, t = .5}))

[0.+0.1334270731e-62*I, 0., 0., 0.3686785e-3]

plot(eval(sol6, t = 0.1e-2), r = 1 .. 2)

Similar to Tom Leslie's numerical solution.

plot3d(sol6, r = 1 .. 2, t = 0 .. 1)



Download spherical3.mw

@Daniel Skoog wrote a package to access OEIS - see his MaplePrimes  post about it.

This works fine in Maple 2017.

Download integr.mw


3 4 5 6 7 8 9 Last Page 5 of 35