dharr

Dr. David Harrington

6611 Reputation

21 Badges

20 years, 73 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 answers submitted by dharr

"I tried converting an outupt to a string which removed all the backslashes in the file path. Maybe there is a dedicated command that I overlooked."

SplitPath does this, in a way that doesn't care whether the separator is \\ or /

Other related manipulations are in the FileTools package

f:=currentdir(); #some path

"C:\Users\dharr\OneDrive - University of Victoria\Desktop"

FileTools:-SplitPath(f);

["C:", "Users", "dharr", "OneDrive - University of Victoria", "Desktop"]

NULL

Download SpliPath.mw

numnew := nops(new): numpast := nops(past):
DataFrame(Matrix(numnew, numpast, (i,j) -> {new[i][]} intersect {past[j][]}),
'rows' = [cat("new ", 1..numnew)], 'columns' = [cat("past ", 1..numpast)]);

Download table.mw

Here's a spruced-up version:

Download table2.mw

restart;

with(DEtools):local gamma;

gamma

verg:= (-delta*eta^2 + alpha*eta)*diff(diff(U(xi), xi), xi) - U(xi)*(2*eta*gamma*theta*(delta*eta - alpha)*U(xi)^2 + eta^2*delta*k^2 + (-alpha*k^2 - 2*delta*k)*eta + 2*k*alpha + delta);

(-delta*eta^2+alpha*eta)*(diff(diff(U(xi), xi), xi))-U(xi)*(2*eta*gamma*theta*(delta*eta-alpha)*U(xi)^2+eta^2*delta*k^2+(-alpha*k^2-2*delta*k)*eta+2*k*alpha+delta)

verg2:=collect(verg/(-delta*eta^2 + alpha*eta),U(xi));

 

-2*eta*gamma*theta*(delta*eta-alpha)*U(xi)^3/(-delta*eta^2+alpha*eta)+(-eta^2*delta*k^2-(-alpha*k^2-2*delta*k)*eta-2*k*alpha-delta)*U(xi)/(-delta*eta^2+alpha*eta)+diff(diff(U(xi), xi), xi)

dsolve(verg2);

U(xi) = c__2*((delta*eta^2*k^2-alpha*eta*k^2-2*delta*eta*k+2*alpha*k+delta)/(-c__2^2*delta*eta^2*gamma*theta+c__2^2*alpha*eta*gamma*theta+delta*eta^2*gamma*theta+eta^2*delta*k^2-alpha*eta*gamma*theta-alpha*eta*k^2-2*delta*eta*k+2*k*alpha+delta))^(1/2)*JacobiSN(((eta*(-delta*eta+alpha)*(-delta*eta^2*gamma*theta-delta*eta^2*k^2+alpha*eta*gamma*theta+alpha*eta*k^2+2*delta*eta*k-2*alpha*k-delta))^(1/2)*xi/(eta*(-delta*eta+alpha))+c__1)*((delta*eta^2*k^2-alpha*eta*k^2-2*delta*eta*k+2*alpha*k+delta)/(-c__2^2*delta*eta^2*gamma*theta+c__2^2*alpha*eta*gamma*theta+delta*eta^2*gamma*theta+eta^2*delta*k^2-alpha*eta*gamma*theta-alpha*eta*k^2-2*delta*eta*k+2*k*alpha+delta))^(1/2), c__2*(-(-delta*eta^2*gamma*theta-delta*eta^2*k^2+alpha*eta*gamma*theta+alpha*eta*k^2+2*delta*eta*k-2*alpha*k-delta)*eta*gamma*theta*(-delta*eta+alpha))^(1/2)/(-delta*eta^2*gamma*theta-delta*eta^2*k^2+alpha*eta*gamma*theta+alpha*eta*k^2+2*delta*eta*k-2*alpha*k-delta))

Notice the JacobiSN, so look this up in DLMF

Imported from https://dlmf.nist.gov/22.13.E13 This is presumably what you mean by a standard form. (Paste the pMML into Maple and change it to 2D math.)

diff(sn(z, k), z, z) = -(k^2+1)*sn(z, k)+2*k^2*sn(z, k)^3

diff(diff(sn(z, k), z), z) = -(k^2+1)*sn(z, k)+2*k^2*sn(z, k)^3

This is close the above, so Maple must have found some change of variables

map(factor,[dcoeffs(verg2,U(xi))]);

[1, 2*gamma*theta, (-delta*eta^2*k^2+alpha*eta*k^2+2*delta*eta*k-2*alpha*k-delta)/(eta*(-delta*eta+alpha))]

NULL

Download ode.mw

Something like this? I generated Cases automatically, but am not sure exactly about the rules.
[Edit: now handles the case where neither A nor B are zero,]

restart

with(PDEtools); with(DEtools)

undeclare(prime)

`There is no more prime differentiation variable; all derivatives will be displayed as indexed functions`

c := combinat:-powerset({A, B}); Cases := [seq({seq(x = 0, `in`(x, s))}, `in`(s, c))]

{{}, {A}, {B}, {A, B}}

[{}, {A = 0}, {B = 0}, {A = 0, B = 0}]

S := diff(G(xi), xi) = G(xi)^2+A*G(xi)+B

diff(G(xi), xi) = G(xi)^2+A*G(xi)+B

"sol:=table():  for i,case in Cases do    sol[i]:=dsolve(eval(S,case));  end do;  "

G(xi) = -(1/2)*A-(1/2)*tanh((1/2)*(A^2-4*B)^(1/2)*(c__1+xi))*(A^2-4*B)^(1/2)

G(xi) = tan(B^(1/2)*(c__1+xi))*B^(1/2)

G(xi) = A/(-1+exp(-A*xi)*c__1*A)

G(xi) = 1/(-xi+c__1)

NULL

Download find_all_case_solution_of_ode_.mw

 

That is strange. Looks like a bug. You can find the singular vector by solving linear equations. (Matrices don't show correctly on Mapleprimes for Maple 2024.)

with(LinearAlgebra); assume(k > 0); assume(abs(l) < 1)

M := Matrix(2, 2, {(1, 1) = l-(1-l*conjugate(l))*k, (1, 2) = -a*(1-l*conjugate(l))*k, (2, 1) = 0, (2, 2) = l+(1-l*conjugate(l))*k})

Matrix(%id = 36893490269467697748)

simplify(SingularValues(M)[1])

(1/2)*(2*(-abs(l)^2+1)*k*(k^2*(abs(l)-1)^2*(abs(l)+1)^2*abs(a)^4+(4*k^2*abs(l)^4+(-8*k^2+4)*abs(l)^2+4*k^2)*abs(a)^2+4*l^2+8*abs(l)^2+4*conjugate(l)^2)^(1/2)+2*k^2*(abs(a)^2+2)*abs(l)^4+4*(-abs(a)^2*k^2-2*k^2+1)*abs(l)^2+2*k^2*(abs(a)^2+2))^(1/2)

X := LinearAlgebra:-HermitianTranspose(M).M; Y := M.LinearAlgebra:-HermitianTranspose(M)

Matrix(2, 2, {(1, 1) = (l-(1-l*conjugate(l))*k)*conjugate(l-(1-l*conjugate(l))*k), (1, 2) = -conjugate(l-(1-l*conjugate(l))*k)*a*(1-l*conjugate(l))*k, (2, 1) = -k*conjugate(a*(1-l*conjugate(l)))*(l-(1-l*conjugate(l))*k), (2, 2) = a*(1-l*conjugate(l))*k^2*conjugate(a*(1-l*conjugate(l)))+(l+(1-l*conjugate(l))*k)*conjugate(l+(1-l*conjugate(l))*k)})

Matrix(%id = 36893490269631762788)

Eigenvectors are zero, but X is not defective (distinct eigenvals)!! Same happens for Y

evalsX, evecsX := Eigenvectors(X); simplify(evalsX[1]-evalsX[2])

Vector[column](%id = 36893490269492440468), Matrix(%id = 36893490269492440588)

(-abs(l)^2+1)*k*(k^2*(abs(l)-1)^2*(abs(l)+1)^2*abs(a)^4+(4*k^2*abs(l)^4+(-8*k^2+4)*abs(l)^2+4*k^2)*abs(a)^2+4*l^2+8*abs(l)^2+4*conjugate(l)^2)^(1/2)

So do it the hard way - we get some arbitrary scale factor this way, which we set to 1

uX1 := simplify(LinearSolve(X-evalsX[1]*IdentityMatrix(2), `<,>`(0, 0), free = t)); uX1 := eval(uX1, t[1] = 1); uX2 := simplify(LinearSolve(X-evalsX[2]*IdentityMatrix(2), `<,>`(0, 0), free = t)); uX2 := eval(uX2, t[1] = 1); uX := `<|>`(uX1, uX2)

Vector(2, {(1) = 1, (2) = (1/2)*((-abs(l)^2+1)*k*sqrt(k^2*(abs(l)-1)^2*(abs(l)+1)^2*abs(a)^4+(4*k^2*abs(l)^4+(-8*k^2+4)*abs(l)^2+4*k^2)*abs(a)^2+4*l^2+8*abs(l)^2+4*conjugate(l)^2)+abs(a)^2*abs(l)^4*k^2+(-2*abs(a)^2*k^2-2*k*l-2*conjugate(l)*k+2)*abs(l)^2+(2*k-2*l)*conjugate(l)+k*(abs(a)^2*k+2*l))/(a*k*(abs(l)-1)*(abs(l)+1)*(k*abs(l)^2-k+conjugate(l)))})

Vector[column](%id = 36893490269491420204)

Check - -2*Re(l)+conjugate(l)+lshould simplify to zero but doesn't.

Xtest := simplify(X.uX-uX.DiagonalMatrix(evalsX))

Matrix(%id = 36893490269566840092)

So try with some random values;

eval(Xtest, {a = 2.1, k = 1.9, l = .3+.2*I})

Matrix(%id = 36893490269505921140)

NULL

Download singular_vect_2_by_2_matrix.mw

See the help page ?encrypted. This means that for procedures in a repository with option encrypted, the code cannot be seen by the user after the procedure is loaded.

Maple has a web page for third party products; not sure how many are packages. As for licensing, I guess you need to read your licencing agreement carefully to see what is allowed.

See also here.

I got a different error (Maple 2024.1). If I take the numerator it is a cubic and can be solved, but the answer is probably not useful, unless you put some numbers in.

Q_simplify_solv.mw

The signs of the terms in the Robin bcs have to be consistent to get a physically meaningful solution. I'm more used to thinkng about diffusion with chemical reactions at the boundary, where the rate of production in the reaction at x=0 has to match the flux in solution in sign and magnitude, which implies that you need a negative sign in your alpha*T(0,t) term if your beta coefficient is to be positive. For your other boundary your temperature is negative (?!) and so that term is already of opposite sign. So here's how to fix it.

Since my problem involving a 2d heat equation was answered so well,

(hats off to nm) I thought I would try another PDE question.

 

Question: can pdsolve handle a combination of the function and its derivative,

i.e., Robin boundary conditions?

 

I looked through the document example, pdsolve_boundaryconditions, which was

beautifully written. However, I saw no example with Robin BCs.

 

"restart; interface(imaginaryunit = i):      StartTemp(x):= `T__0`*(e)^(-`alpha__0` )*cos(Pi*(x)/(`L__0`)):   `L__0` := 10: #` Length of rod` `T__0` := 100: #` max temperature` `alpha__0` := 6/(100): #`  Dampening factor` plot(StartTemp(x) , x = 0 .. `L__0`,  thickness = 5,             view = [0..`L__0`, -100..100],  labels = ["Position on rod", "Temperature"], labelfont = [Times, 14],             labeldirections = [horizontal, vertical], title = "Temperature in rod at t = 0" , titlefont = [Times, 16],             size = [500, 225]);"

The equation and the initial condition:

heat_equation := diff(T(x, t), t) = k*(diff(T(x, t), x, x)); initial_condition := T(x, 0) = T__i(x)

Using a linear combination of the function and its derivative, i.e., Robin BCs. For positive beta, need first term negative at x = 0

boundary_conditions := -`&alpha;__r`*T(0, t)+`&beta;__r`*(D[1](T))(0, t) = f(t), `&alpha;__r`*T(L__0, t)+`&beta;__r`*(D[1](T))(L__0, t) = g(t)

 

Setting the constants

"  `alpha__r` := 1:  #` Constant for boundary condition`  `beta__r` := 2: #` Constant for boundary condition`     `T__s`:= 10:  #` Time interval`  k:=9/(10):  #` Heat-conductivity of material`    f(t) :=-1/(10) t:                  g(t):=1/(10) t:  "

Solving using nm's technique of not defining the function before calling pdsolve

We now have an eigenvalue type solution, for which an analytical solution is not available for the eigenvalues.

sols1 := pdsolve({heat_equation, boundary_conditions, initial_condition}, T(x, t))

T(x, t) = `casesplit/ans`(-10/9+(1/18)*x^2-(5/9)*x+Sum((lambda[n]*cos((1/10)*lambda[n]*x)+5*sin((1/10)*lambda[n]*x))*exp(-(9/1000)*lambda[n]^2*t)*lambda[n]*(Int(-(x^2-18*T__i(x)-10*x-20)*(lambda[n]*cos((1/10)*lambda[n]*x)+5*sin((1/10)*lambda[n]*x)), x = 0 .. 10))/((45*lambda[n]^2-1125)*sin(2*lambda[n])+90*lambda[n]*(lambda[n]^2-5*cos(2*lambda[n])+30)), n = 0 .. infinity)+(1/10)*t, {And(tan(lambda[n])*lambda[n]^2-25*tan(lambda[n])-10*lambda[n] = 0, 0 < lambda[n])})

Now include the initial temperature function

sols := eval(sols1, T__i(x) = StartTemp(x))

Extract the pieces we need. (value works out the integral, which fortunately has an analytical solution.)

summand, rest := selectremove(has, op(1, rhs(sols)), Sum); summand := value(op(1, summand)); eigenvalueeqn := lhs(op(1, op(2, rhs(sols))[])); eigenvalname := indets(eigenvalueeqn, name)[]; eigenvalfn := MakeFunction(eigenvalueeqn, eigenvalname)

200*(lambda[n]*cos((1/10)*lambda[n]*x)+5*sin((1/10)*lambda[n]*x))*exp(-(9/1000)*lambda[n]^2*t)*(90*exp(-3/50)*lambda[n]^5*sin(lambda[n])+Pi^2*sin(lambda[n])*lambda[n]^3-450*exp(-3/50)*cos(lambda[n])*lambda[n]^4-sin(lambda[n])*lambda[n]^5-10*Pi^2*cos(lambda[n])*lambda[n]^2-450*exp(-3/50)*lambda[n]^4+10*cos(lambda[n])*lambda[n]^4-15*Pi^2*sin(lambda[n])*lambda[n]+15*sin(lambda[n])*lambda[n]^3-50*Pi^2*cos(lambda[n])+50*lambda[n]^2*cos(lambda[n])+50*Pi^2-50*lambda[n]^2)/(lambda[n]^2*(Pi^2-lambda[n]^2)*((45*lambda[n]^2-1125)*sin(2*lambda[n])+90*lambda[n]*(lambda[n]^2-5*cos(2*lambda[n])+30)))

tan(lambda[n])*lambda[n]^2-25*tan(lambda[n])-10*lambda[n]

lambda[n]

proc (y1) options operator, arrow; tan(y1)*y1^2-25*tan(y1)-10*y1 end proc

plot(eigenvalfn, 0 .. 20)

So we need to numerically find the eigenvalues and assemble the solution. Using NextZero with default settings is unreliable. Sturm-Liouville theory can be used to bracket the roots (see here for an example), but by inspection it is clear that the roots lie between successive multiples of Pi, which helps fsolve find them reliably.

nt := 100; # number of terms
r := fsolve(eigenvalfn(x), x = 0.1..3):
sol1 := eval(summand, eigenvalname = r):
for j from 2 to nt do
  r := fsolve(eigenvalfn(x),x=(j-1)*Pi..j*Pi);
  sol1 := sol1 + eval(summand, eigenvalname = r);
  #print(j, r, evalf(j*Pi));
end do:
T__sol := MakeFunction(sol1 + rest, x, t):

100

In plotting the solution at t = 0, we now reproduce the initial temperature function.

plot(T__sol(x, 0), x = 0 .. L__0, -120 .. 120, thickness = 5, size = [500, 225])

 

NULL

Download Heat_Robin.mw

Do you want b to be an integer? Maple's RootOf for polynomials (where all the powers are integers can be quite useful), but for more complicated cases like here where b is not necessarily integer, there are probably only a few things one can do. The attached explores some of the parameter space.

Roots_of_a_Polynomial.mw

When you solved with the inequalities, solve returned a RootOf with a range to say which of the roots to use, which was the one evalf produced, and which presumably satisfies the equations and inequalities.

When you solved with just the equations, solve returned a generic RootOf, which stands for all roots, but evalf gives only the first. If you use allvalues and then evalf, you will get all possibilities, and the 7th one is the one that agrees with the earlier solve.

problems_with_solve_15.10.24.mw

iroot is much faster, but it gives the closest integer, not the floor, so there may be a digit difference.

Edit: floor(root(evalf(n), 3)) is faster yet, but needs appropriately high Digits.

restart;

Digits:=50;

50

n := 33333333333333333333333333333333333444444444444444444;
CodeTools:-Usage(floor(root(evalf(n), 3)),iterations=1000);
CodeTools:-Usage(floor(root(n, 3)),iterations=1000);
CodeTools:-Usage(iroot(n, 3, 'exact'),iterations=1000);
exact;

33333333333333333333333333333333333444444444444444444

memory used=1.75KiB, alloc change=0 bytes, cpu time=16.00us, real time=16.00us, gc time=0ns

321829794868543252

memory used=29.20KiB, alloc change=0 bytes, cpu time=172.00us, real time=355.00us, gc time=0ns

321829794868543252

memory used=7.45KiB, alloc change=28.99MiB, cpu time=47.00us, real time=85.00us, gc time=46.88us

321829794868543253

false

Check

Digits := 50;
root(evalf(n),3);

50

321829794868543252.61997759481168897289518515510196

 

Download iroot.mw

Since you want your answer in terms of h1 and h2, you want to solve for the other variables.

ans := solve({v1^2 = 2*g*(h-h1), (1/2)*g*t^2 = h2, v1*t+(1/2)*g*t^2 = h1}, {g, h, v1})

{g = 2*h2/t^2, h = (1/4)*(h1^2+2*h1*h2+h2^2)/h2, v1 = (h1-h2)/t}

simplify(ans)

{g = 2*h2/t^2, h = (1/4)*(h1+h2)^2/h2, v1 = (h1-h2)/t}

NULL

Download Expression_of_NLP_with_desired_parameters.mw

If you substitute the solve solutions back into the equation, you find that (except for the trivial solution), they are not actually solutions. On the other hand, fsolve gives an accurate solution for one range I tried (I was assuming you wanted a positive solution). Since you want a numerical solution, fsolve is preferred.

By the way, your second equation is c__3*(expression), so c__3 = 0 always works. If that is a solution you are not interested in, then I would change that to just (expression); then it may be easier to get a solution without c__3 = 0.

question11.mw

Aside from the 2I, there is a 6I (interpreted correctly anyway), and then missing spaces after a C1 and an h3 which makes them into unknown functions. Then since everything else is real you can use evalc() or evalc(Re())

real_and_imaginary_.mw

As the others have pointed out the sum involves roots of a degree 4 polynomial. If you put some values in, then this will give you directly what you want. (I used some random values so didn't get realistic output.)

Download RootOf.mw

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