dharr

Dr. David Harrington

6728 Reputation

21 Badges

20 years, 98 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

@janhardo The DLMF (and the corresponding hardback book version) is produced by NIST (National Institute of Standards and Technology) in the U,S., and is the successor to Abramowitz and Stegun's, "Handbook of Mathematical Functions"; it is considered an authoritative source, perhaps "the" authoritative source.

These differential equations are standard forms in the sense that they are often quoted, e.g., one speaks of "Bessel's equation"  and there is a common way to write it, though of course there is not universal consistency in these things. 

But in terms of classifications in a text as I think you want, there is probably not much of a standard. I personally like Zwillinger's "Handbook of Differential Equations". Of course these classifications exist because of different solving methods, but since the objective is often a special function, they are related to the other standard forms.

Maple' odeadvisor(verg) gives: [[_2nd_order, _missing_x], [_2nd_order, _reducible, _mu_x_y1]]

Since many of the second-order differential equations and their "special" functions are related through transformations to the Sturm-Liouville eqation, perhaps it is the standard form to rule them all :-).

@MichaelVio Here's an example.

restart

with(PDEtools):

Differential equation (=0)

det := diff(E(t),t) - A*E(t)/(exp(E(t)/k/T)-1);

diff(E(t), t)-A*E(t)/(exp(E(t)/(k*T))-1)

dchange({t=1/x},det);

-(diff(E(x), x))*x^2-A*E(x)/(exp(E(x)/(k*T))-1)

dchange({t=1/x},det/x^2,expand);

-(diff(E(x), x))-A*E(x)/(x^2*(exp(E(x)/(k*T))-1))

NULL

Download dchange.mw

@salim-barzani I'm not sure I fully understand the order you want to do things or what you want to store, but perhaps this:

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

Fode := (-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)

case1 := {alpha = delta*(A^2*eta^2-2*eta^2*k^2-4*B*eta^2+4*eta*k-2)/(A^2*eta-2*eta*k^2-4*B*eta+4*k), eta = eta, a[0] = -A/(2*gamma*theta*RootOf(_Z^2*gamma*theta+1)), a[1] = RootOf(_Z^2*gamma*theta+1)}

{alpha = delta*(A^2*eta^2-2*eta^2*k^2-4*B*eta^2+4*eta*k-2)/(A^2*eta-2*eta*k^2-4*B*eta+4*k), eta = eta, a[0] = -(1/2)*A/(gamma*theta*RootOf(_Z^2*gamma*theta+1)), a[1] = RootOf(_Z^2*gamma*theta+1)}

K := U(xi) = a[0]+a[1]*G(xi)

U(xi) = a[0]+a[1]*G(xi)

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

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

"for i,case in Cases do    sol:=dsolve(eval(S,case));    case1AB:=eval(case1,case):     Usol:=eval(eval(K,{sol} union case1AB),case);   simplify(odetest(Usol,eval(Fode,case1AB)));  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)

{alpha = delta*(A^2*eta^2-2*eta^2*k^2-4*B*eta^2+4*eta*k-2)/(A^2*eta-2*eta*k^2-4*B*eta+4*k), eta = eta, a[0] = -(1/2)*A/(gamma*theta*RootOf(gamma*theta*_Z^2+1)), a[1] = RootOf(gamma*theta*_Z^2+1)}

U(xi) = -(1/2)*A/(gamma*theta*RootOf(gamma*theta*_Z^2+1))+RootOf(gamma*theta*_Z^2+1)*(-(1/2)*A-(1/2)*tanh((1/2)*(A^2-4*B)^(1/2)*(c__1+xi))*(A^2-4*B)^(1/2))

0

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

{alpha = delta*(-2*eta^2*k^2-4*B*eta^2+4*eta*k-2)/(-2*eta*k^2-4*B*eta+4*k), eta = eta, a[0] = 0, a[1] = RootOf(gamma*theta*_Z^2+1)}

U(xi) = RootOf(gamma*theta*_Z^2+1)*tan(B^(1/2)*(c__1+xi))*B^(1/2)

0

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

{alpha = delta*(A^2*eta^2-2*eta^2*k^2+4*eta*k-2)/(A^2*eta-2*eta*k^2+4*k), eta = eta, a[0] = -(1/2)*A/(gamma*theta*RootOf(gamma*theta*_Z^2+1)), a[1] = RootOf(gamma*theta*_Z^2+1)}

U(xi) = -(1/2)*A/(gamma*theta*RootOf(gamma*theta*_Z^2+1))+RootOf(gamma*theta*_Z^2+1)*A/(-1+exp(-A*xi)*c__1*A)

0

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

{alpha = delta*(-2*eta^2*k^2+4*eta*k-2)/(-2*eta*k^2+4*k), eta = eta, a[0] = 0, a[1] = RootOf(gamma*theta*_Z^2+1)}

U(xi) = RootOf(gamma*theta*_Z^2+1)/(c__1-xi)

0

NULL

NULL

Download find_all_case_solution_of_ode_.mw

For Vectors, use a:=<2,1>, not <<2,1>>, which is a 2x1 Matrix. Then you can use "." for the dot product.

@Carl Love I think shape = symmetrix or shape = hermitian change the algorithm for float matrices to produce only real eigenvalues and eigenvectors but have no effect in the symbolic case.

@rzarouf You're welcome. You found a Maple bug, and I have submitted an SCR (bug report), so it is fixed in future versions.

@Carl Love The method is fine, but M^%H.M is not symmetric, so using shape=symmetric forces the offdiagonals equal. So you should omit that.

@Scot Gould Although @Rouben Rostamian has clarified that the any Robin boundaries can be solved, it seems that Maple finds sometimes the right answer consistent for the initial condition given and sometimes not. Your first solution was already compatible for the right-hand boundary (at 50 points it is pretty close on almost all the right half of the domain). So I knew to play around with the left-hand one, and used an analogy with the solutions that I was most familar with to figure out that sign change. I usually solve cases that evolve toward steady state, so I added my bias there. But maybe Maple is also consistently solving some classes of boundary conditions; it would help if we know more clearly what assumptions it is making. In this respect, @nm's suggestion to put in the initial solution in afterwards might be counterproductive; maybe Maple is not returning a solution because it knows it can't solve the case compatible with that initial conditions. I suspect very specific assumptions on signs of parameters at the pdsolve stage might be helpful.

I don't think the explanation on the Wikipedia page is that helpful; one needs to think of the relative signs to achieve the right physical outcome. I deliberately kept the fluxes the same sign to get to a steady-state, then adjusted the other signs. The general steady state solution is a straight line. If the slope of that is positive, the flux is negative and there will be the same heat flux coming in the right-hand side as going out the left-hand side at steady state. So if the initial fluxes are not equal but of the correct sign, then the other parts of the Robin boundary condition have to allow the slopes to evolve to be the same (assuming you want steady-state of course). In your case the fluxes evolve to be zero.

I didn't really want to use all that opcrap to extract the pieces and did try to do it without a loop first, but I also didn't want to spend to much time optimizing it. Make a procedure that determines whether it is an eigenvalue solution without an analytical solution for the eigenvalues might be worthwhile, but in that case I would probably solve it "by hand". Actually with Maple that is not too difficult. If you call pdsolve with HINT=f(t)*g(x), then extracting the ode for t and x and relating them to the eigenvalue is straightfoward. The second-order ode in x can be solved for {g(x),lambda} and if there is an analytical formula for lambda, Maple will find it; otherwise you have to find the eigenvalues numerically.

I should say, I almost never solve the heat/diffusion equation by this eigenvalue method, because the error functions and x/sqrt(k*t) usually found don't easily come out of it. Instead I find the solution by Laplace transforms. (To get the general solution containing x/sqrt(k*t) you need to use SimilaritySolutions.)

@Rouben Rostamian  Thanks for the clarification. (Yes, I generally expect solutions to evolve toward a steady state.)

@Teep "Do you think that it is possible to extract a general expression for the roots of f (or g) in terms of a and b?" No, if Maple could find that it would not have returned the RootOf, which is the "next best thing" when it can't find an explicit solution. Sometimes it is possible to manually rearrange things to help Maple out, but I think with a non-integer b in the exponent, it is likely an intractable problem. Even for the integer b case, polynomials of high degree do not usually have explicit symbolic solutions.

@Saalehorizontale As you expect, a RootOf(..., 0.35..0.36) or RootOf(...,index=2) corresponds to a single root. Unfortunately, as you already saw, RootOf without an index or range or label sometimes means any/all of the roots and sometimes means the first one. Solving a polynomial gives a RootOf that says all possible roots work, and can be found by allvalues, but evalf just produces one. Other commands, particularly those where a RootOf stands for an algebraic number in a field mean also the first one, e.g. evala(Minpoly(..)).

I didn't check the timings, but it's not unexpected that figuring out the actual numerical values of the roots is slow. On the other hand solving for polynomial roots is optimized compared to roots of more general equation, e.g., fsolve(polynomial in x,x,complex) will quickly return all roots.

@dharr I edited the above to include floor(root(evalf(n), 3)). In floor(root(n, 3)), root(n,3) just gives n^(1/3) and floor does all the work.

@GFY You can solve it by narrowing the ranges, with guidance of a plot. You could modify the equations so a__3=0,c__3=0 is not a solution by factoring and dividing through by power(s) of a__3 and c__3.

question11.mw

@GFY That was unexpected for me. Solve is really for symbolic solutions, and you have floats in your equations. I always formulate equations symbolically, and try solve. Then one can substitute values for the parameters, e.g., params:= {a=3.1, b=2.1 etc}, then eval(eqn, params);

@Aixleft math Typically I specify as many variables as equations, but if there is a solution with more, Maple can often find it.

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