Scot Gould

Scot Gould

692 Reputation

14 Badges

11 years, 27 days
Claremont McKenna, Pitzer, Scripps College
Professor of Physics
Upland, California, United States
Dr. Scot Gould is a professor of physics in the W.M. Keck Science Department of Claremont McKenna, Pitzer, and Scripps Colleges - members of The Claremont Colleges in California. He was involved in the early development of the atomic force microscope. His research has included numerous studies and experiments using scanning probe microscopes, particularly those involving natural fibers such as spider silk. More recently, he was involved in developing and sustaining AISS. This full-year multi-unit, non-traditional, interdisciplinary undergraduate science education course integrated topics from biology, chemistry, physics, mathematics, and computer science. His current interest is integrating computational topics into the physics curriculum. He teaches the use of Maple's computer algebraic and numerical systems to assist students in modeling and visualizing physical and biological systems. His Dirac-notation-based quantum mechanics course is taught solely through Maple.

MaplePrimes Activity

These are answers submitted by Scot Gould

I see two problems.


1) `D` is reserved in Maple as an operator. So, use a different constant, such as k,

for this heat-equation type PDE.


2) Are you sure about those conditions? If I ask pdsolve to provide a solution

without any conditions, it returns two ODEs, one of which requires an initial

condition and the other requires only two boundary conditions.


pdsolve(diff(u(z, t), t) = k*(diff(u(z, t), z, z))+A*z*exp(lambda*z))

PDESolStruc(u(z, t) = f__1(z)*f__2(t)-((lambda*z-2)*c__1*exp(lambda*z)/lambda^3+c__2*z+c__3)*A/(k*c__1), [{diff(f__2(t), t) = k*_c[1]*f__2(t), diff(diff(f__1(z), z), z) = _c[1]*f__1(z)}])





Calculating the coefficients as given by the equation for a generic function f of

the variable x. (BesselJ(0, x) is created by writing BesselJ when needed, hitting the esc button, selecting

the proper character, and filling in the details.)


"A(n, f, x) := ((∫)[0]^(1.0)x*f*(J)[0](alpha(n)*x) ⅆx)/((∫)[0]^(1.0)x*(J)[0](alpha(n)*x)^(2) ⅆx) :"



"  alpha(n) := BesselJZeros(0, n): "


Testing for an example, such as yours.


"`f__op`(x) := (1-x^(2)):    '`A__3`' = A(3,`f__op`(xi), xi)"

A__3 = 0.4547647069e-1



Here is an example of how adding more terms does a better job in the

approximation. The function approximates up to N terms.

" `f__Bessel`(f, x, N) := (∑)( A(n, f, x)*(J)[0](alpha(n)*x))  :"


Slide the slider to show how many terms are used to approximate:


Explore(plot([f__op(xi), f__Bessel(f__op(xi), xi, N)], xi = 0 .. 1, 0 .. 1, legend = [f, f__Bessel], size = [400, 400]), parameters = [N = 1 .. 5], initialvalues = [N = 2])




Here is a worksheet I give new users that provides a simple example of animating ODEs using the Explore procedure.

@dharr 's answer looks fine to me



Let me give you one other example that might help:


A is a random matrix


A := LinearAlgebra:-RandomMatrix(6, 7, generator = 0 .. 10)

Matrix(%id = 36893489902322627876)


M is a matrix of 999s


M := Matrix(3, 4, 999)

Matrix(%id = 36893489902322619924)



Into A, from rows 3 to row 5 at column 2 through column 3, copy M,

rows 1 through rows 3 at columns 3 through column 4.


A[3 .. 5, 2 .. 3] := M[1 .. 3, 3 .. 4]

Matrix(%id = 36893489902322617148)


Now, what is A?


Matrix(%id = 36893489902322627876)






I'm a bit confused as to what you want to do. Extend is a procedure in the ArrayTools package.


Are you trying to add another column?  If so, two elements are required. Are you trying to add

another row? That requires four elements.  But you provided three. Hence the confusion.

Personally, I prefer to use the brackets technique to create vectors or matrices:


restart; XX := `<,>`(`<|>`(1, 1, 1, 1), `<|>`(2, 2, 2, 2))

Matrix(%id = 36893490419552086124)



Add another column of zeros by adding a vector `<,>`(0, 0) to the matrix


XX2 := `<|>`(XX, `<,>`(0, 0))

Matrix(%id = 36893490419552074796)


Add a ";" means add another row.

XX2 := `<,>`(`<|>`(XX), `<|>`(0, 0, 0, 0))

Matrix(%id = 36893490419552063844)




Just a heads up. There is a problem in your worksheet that others did not point out,

although mmcdara corrected.


When you extract your solution, you cannot immediately apply it as the expression

for a function. (I wish one could, but it would break Maple.)   It may look like it is

working, but it is not.


Referring to your example, and after adding in the value procedure: (Allow me to

use 2d-input. This input does not affect the point of this post.) 


restart; N := 4; T := 5.0*(1/60); M := 6905

"dose(t) := M*(&sum;)(Heaviside(t-(24*k)/(N))*(e)^((((-t+(24*k)/(N)))/(T)))) :  "

deq := diff(C(t), t) = dose(t)-Cl*C(t)


Going through step-by-step


deq_sol := dsolve({deq, C(0) = 0}); Cl := .32


Simplifying deq_sol using the value procedure:


deq_sol := `assuming`([value(deq_sol)], [Cl > 0])


Now the problem: Here is what happens when you directly assign the

the expression on the rhs of the equation to a function name:


sol := proc (t) options operator, arrow; rhs(deq_sol) end proc

proc (t) options operator, arrow; rhs(deq_sol) end proc


Let's look at a value for your solution such as one at
t = 3.2




Notice that the value of 3.2 is not substituted into the function.  Instead, Maple

returns the expression. Because it is just the expression, 

the plot procedure works even though this is not a true function.


To make a function, you can either use the older procedure, unapply, or, since

you are using Maple 2024, you can use the new name: MakeFunction.


sol := MakeFunction(rhs(deq_sol), t)

proc (t) options operator, arrow; 6905*(-0.8561643836e-1*Heaviside(t)*exp(-11.68000000*t)+0.8561643836e-1*Heaviside(t)-0.8561643836e-1*Heaviside(t-42.)*exp(504.-11.68000000*t)+58813.22494*Heaviside(t-42.)-0.8561643836e-1*Heaviside(t-36.)*exp(432.-11.68000000*t)+8622.428242*Heaviside(t-36.)-0.8561643836e-1*Heaviside(t-30.)*exp(360.-11.68000000*t)+1264.108011*Heaviside(t-30.)-0.8561643836e-1*Heaviside(t-24.)*exp(288.-11.68000000*t)+185.3270353*Heaviside(t-24.)-0.8561643836e-1*Heaviside(t-18.)*exp(216.-11.68000000*t)+27.17023364*Heaviside(t-18.)-0.8561643836e-1*Heaviside(t-12.)*exp(144.-11.68000000*t)+3.983345414*Heaviside(t-12.)-0.8561643836e-1*Heaviside(t-6.)*exp(72.-11.68000000*t)+.5839861703*Heaviside(t-6.))*exp(-.32*t) end proc


Now it will act like a function:





Area_under_curve := int(sol(t), t = 0 .. 48.0)




I also agree with mmcdara that you would benefit from creating a function out of both

t and Cl.


    sol := MakeFunction(rhs(deq_sol), t, Cl)


I hope you found my comment useful.




I concur. It appears the only way to generate a solution is to do

it numerically. However, one can make a function out of the integral.


If we pull out the term: f^s/factorial(s) , we are left with just the integral

"S(n,s) := (&int;)[0]^(1.0)((n*p+(e)^(n*p))*(e)^(-p)*p^(s-1))/((n*p+(e)^(2*n*p)) ) &DifferentialD;p :"


Testing it for n = 2, s = 3:

S(2, 3)*f^3/factorial(3)



I don't have time to try it out, but you might be able to generate

a polynomial series approximation to the integral from the data.

Looking at two plots.


plot(S(n, 3), n = 0 .. 5, 0 .. S(0, 3), size = [500, 200])


plot3d(S(n, s), n = 0 .. 5, s = 1 .. 5, labels = [n, s, S], labelfont = [Times, 20], size = [500, 500], orientation = [64, 74, -2])




Since you're trying to make a function, use the procedure MakeFunction. It works the same way as unapply and is more intuitive and memorable.

@Ronan, I'm not sure what you are asking about sorting based on the index of min value, but as nm says, and what I do for the reason you are experiencing, I do sort and maintain the order of the eigenvalues. The minimum is the first one.

I can recommend this video where Dr. Gerhard of Maplesoft discusses writing efficient code. It includes a few useful comments about the CodeTools package.

They have also provided the document Dr. Gerhard uses in the presentation:


Given its elegance, I do prefer using the elementwise operator (~),

as nicely noted by others. However, if you want something that clearly shows

you are creating a vector, i.e., using the vector brackets < >,  may I recommend

the seq procedure, which is possibly the most useful procedure in Maple :


"h(t,x) := ((m*t^(2)+6*t-2*x())^(2))/(36*g*t^(2)) :  `t__values`:=1, 5, 10, 15, 20:  "

V := `<,>`(seq(h(t, x), t = t__values))

Vector[column](%id = 36893489650694321204)





Looking at the expression:

expression := 5*x*cos((1/2)*x^2)-5*sin(5*x*(1/2))*cos(5*x*(1/2)); plot(expression, x = -5 .. 5)


It is pretty clear that with the cosines and sines, there are an infinite number

of solutions. So, this limits you to finding only "floating point" solutions, i.e.,

only numerical solutions. That is why your choice of "fsolve" makes sense.


Here are the first 4 solutions for positive x

sols := fsolve(expression = 0, x = 0 .. 5, maxsols = 4)

.4258227298, 1.580743250, 3.084753386, 3.939229633


To which you can refer to each of them by an index

solution1 := sols[1]





restart; sols := dsolve({diff(y(x), x) = y(x)-x, diff(z(x), x) = -x-y(x), y(0) = 0, z(0) = 0})

{y(x) = x+1-exp(x), z(x) = -x^2-x+exp(x)-1}


y_expression := eval(y(x), sols)



y := MakeFunction(y_expression, x)

proc (x) options operator, arrow; x+1-exp(x) end proc


z := MakeFunction(eval(z(x), sols), x)

proc (x) options operator, arrow; -x^2-x+exp(x)-1 end proc


plot([y(x), z(x)], x = 0 .. 1)





restart; interface(imaginaryunit = I)

n := 3; m := 2; eqx := x^(n/m) = a; maple_sol := [PDEtools:-Solve(eqx, x)]

[x = a^(2/3), x = (1/4)*a^(2/3)*(1+I*3^(1/2))^2, x = (1/4)*a^(2/3)*(I*3^(1/2)-1)^2]


test_eqns := `~`[eval](eqx, maple_sol)

[a = a, (1/16)*4^(1/2)*(a^(2/3)*(1+I*3^(1/2))^2)^(3/2) = a, (1/16)*4^(1/2)*(a^(2/3)*(I*3^(1/2)-1)^2)^(3/2) = a]


For "a = 1, " only the first solution works:

`~`[is](`~`[eval](test_eqns, a = 1))

[true, false, false]


For a = -1, two solutions work

`~`[is](`~`[eval](test_eqns, a = -1))

[true, true, false]


For a = I, a different pair of solutions work.

`~`[is](`~`[eval](test_eqns, a = -I))

[true, false, true]





I like this problem because it forced me to think about what

happens in the expression x^(1/3).


Since Maple knows nothing about x, it could be complex. Hence,

let's rewrite x asa+I*b, where a and b are real numbers.


(Note: For this worksheet, since I find the letter 'i' more readable

than 'I', there is a change in the interface.)


restart; interface(imaginaryunit = I, showassumed = 0)


And here is where we define x. We tell Maple to assume that a, b 

are always real.


x := a+I*b; assume(a::real, b::real)


Let's start with the easy test. Is (x^(1/3))^3 = x ?

(x^(1/3))^3 = x

a+I*b = a+I*b


Yep, looks good.


Now, we try again with the other equation. Is (x^3)^(1/3) = x ?


eq := (x^3)^(1/3) = x

((a+I*b)^3)^(1/3) = a+I*b



Hmm, not good. Let's see if we can get Maple to expand the equation


eq := expand(eq)

(a^3+(3*I)*a^2*b-3*a*b^2-I*b^3)^(1/3) = a+I*b


It looks like the left-hand side of the equation is not equal to the right-hand side.

Pull out the trusty evalc to see if we can do better. The procedure evalc 

evaluates the expression, assuming that all variables in the expression are real

and then collects the real and imaginary parts separately.


eq := evalc(eq)

(a^2+b^2)^(1/2)*cos((1/3)*arctan(3*a^2*b-b^3, a^3-3*a*b^2))+I*(a^2+b^2)^(1/2)*sin((1/3)*arctan(3*a^2*b-b^3, a^3-3*a*b^2)) = a+I*b


Definitely not the same. The real part of the left-hand side is not equal to a,

nor do the imaginary parts equate. Hence, this is why Maple cannot

and should not simplify (x^3)^(1/3) to x.


But what if we say we know x is real? Then, there is no imaginary

component. Hence, b = 0.


eq := eval(eq, b = 0)

(a^2)^(1/2)*cos((1/6)*(1-signum(0, a, 1))*Pi)+I*(a^2)^(1/2)*sin((1/6)*(1-signum(0, a, 1))*Pi) = a


Still not the same. However, there are two conditions if x is real. Either a is

negative or a is positive.  So, let's simplify for both cases:


`assuming`([simplify(eq)], [a < 0])

-(1/2)*a*(1+I*3^(1/2)) = a


Nope, the equation is not true even if x is real


For a > 0 :

`assuming`([simplify(eq)], [a >= 0])

a = a




Conclusion (x^3)^(1/3) = x if and only if x >= 0"." 

That is, only if x is a nonnegative real number.



1 2 3 4 5 6 Page 1 of 6