ecterrab

14645 Reputation

24 Badges

20 years, 152 days

MaplePrimes Activity


These are answers submitted by ecterrab

@shingy@Carl Love 

Carl is correct. In addition: your equ has r^beta, which makes difficult the differential elimination process (equivalent to a triangularization towards decoupling the PDE system and simplifying it taking into account integrability conditions). So, I didn't wait to see if it returns in reasonable time without using the option integrabilityconditions= false but it returns right-away with this option.

So here are some more details on your problem, including the actual form of the symmetry (not just of the determining PDE system for it) and actual invariant solutions for equ derived directly from these symmetries (although involving an unsolved 2nd order linear ODE that requires specializing beta to get solved to-the-end)

Use this to simplify in significant ways the display of formulas - note that derivatives will then be displayed indexed.

PDEtools:-declare(V(t, r, S), (_xi, _eta)(t, r, S, V), _Y(r))

V(t, r, S)*`will now be displayed as`*V

 

_xi(t, r, S, V)*`will now be displayed as`*_xi

 

_eta(t, r, S, V)*`will now be displayed as`*_eta

 

_Y(r)*`will now be displayed as`*_Y

(1)

equ := diff(V(t, r, S), t)+(1/2)*sigma^2*S^2*(diff(V(t, r, S), S, S))+S*(diff(V(t, r, S), S))-r*V(t, r, S)+(1/2)*Sigma^2*(r^beta)^2*(diff(V(t, r, S), r, r))+(a-`κr`-`λΣ`*r^beta)*(diff(V(t, r, S), r)) = 0

diff(V(t, r, S), t)+(1/2)*sigma^2*S^2*(diff(diff(V(t, r, S), S), S))+S*(diff(V(t, r, S), S))-r*V(t, r, S)+(1/2)*Sigma^2*(r^beta)^2*(diff(diff(V(t, r, S), r), r))+(a-`κr`-`λΣ`*r^beta)*(diff(V(t, r, S), r)) = 0

(2)

The Determining PDE, without using integrability conditions to avoid a difficult differential elimination process involving symbolic powers of r

PDEtools:-DeterminingPDE(equ, integrabilityconditions = false)

{(S*(diff(_xi[r](t, r, S, V), r))-S*(diff(_xi[S](t, r, S, V), S))+_xi[S](t, r, S, V))*r^(1+2*beta)-beta*r^(2*beta)*_xi[r](t, r, S, V)*S, (diff(diff(_xi[t](t, r, S, V), S), V))*S^2*sigma^2-2*(diff(_xi[S](t, r, S, V), V)), -2*r^(1+2*beta)*(diff(_xi[S](t, r, S, V), r))*S*Sigma^2-2*(diff(_xi[r](t, r, S, V), S))*r*S^3*sigma^2, diff(diff(_xi[r](t, r, S, V), V), r)-(1/2)*(diff(diff(_eta[V](t, r, S, V), V), V)), r^(1+2*beta)*(diff(diff(_xi[S](t, r, S, V), V), r))*Sigma^2+(diff(diff(_xi[r](t, r, S, V), S), V))*S^2*r*sigma^2+2*(`λΣ`*r^(1+beta)-r*(a-`κr`))*(diff(_xi[S](t, r, S, V), V)), (diff(diff(_xi[S](t, r, S, V), S), V))*S*sigma^2-(1/2)*(diff(diff(_eta[V](t, r, S, V), V), V))*S*sigma^2-2*(diff(_xi[S](t, r, S, V), V)), r^(1+2*beta)*(diff(diff(_eta[V](t, r, S, V), r), r))*S*Sigma^2+S^3*r*(diff(diff(_eta[V](t, r, S, V), S), S))*sigma^2-2*S*(`λΣ`*r^(1+beta)-r*(a-`κr`))*(diff(_eta[V](t, r, S, V), r))+2*(S*r*V*(diff(_eta[V](t, r, S, V), V))+S^2*(diff(_eta[V](t, r, S, V), S))+S*(diff(_eta[V](t, r, S, V), t))-2*S*(diff(_xi[S](t, r, S, V), S))*r*V-S*r*_eta[V](t, r, S, V)-V*(S*_xi[r](t, r, S, V)-2*_xi[S](t, r, S, V)*r))*r, -(diff(diff(_xi[t](t, r, S, V), r), r))*r^(1+2*beta)*S*Sigma^2-S^3*r*(diff(diff(_xi[t](t, r, S, V), S), S))*sigma^2+2*S*(`λΣ`*r^(1+beta)-r*(a-`κr`))*(diff(_xi[t](t, r, S, V), r))-2*r*(S*(diff(_xi[t](t, r, S, V), V))*V*r+S^2*(diff(_xi[t](t, r, S, V), S))+S*(diff(_xi[t](t, r, S, V), t))-2*S*(diff(_xi[S](t, r, S, V), S))+2*_xi[S](t, r, S, V)), r^(1+2*beta)*(diff(diff(_xi[S](t, r, S, V), r), r))*Sigma^2-2*(diff(diff(_eta[V](t, r, S, V), S), V))*S^2*r*sigma^2+(diff(diff(_xi[S](t, r, S, V), S), S))*S^2*r*sigma^2+(-2*`λΣ`*r^(1+beta)+2*r*(a-`κr`))*(diff(_xi[S](t, r, S, V), r))-2*r*(-3*(diff(_xi[S](t, r, S, V), V))*V*r+S*(diff(_xi[S](t, r, S, V), S))-_xi[S](t, r, S, V)-(diff(_xi[S](t, r, S, V), t))), -S*Sigma^2*(diff(diff(_xi[r](t, r, S, V), r), r)-2*(diff(diff(_eta[V](t, r, S, V), V), r)))*r^(1+2*beta)-r*(diff(diff(_xi[r](t, r, S, V), S), S))*S^3*sigma^2+2*(`λΣ`*r^(1+beta)-r*(a-`κr`))*S*(diff(_xi[r](t, r, S, V), r))-4*(`λΣ`*r^(1+beta)-r*(a-`κr`))*S*(diff(_xi[S](t, r, S, V), S))-2*(diff(_xi[r](t, r, S, V), V))*V*r^2*S-2*(diff(_xi[r](t, r, S, V), S))*r*S^2-2*(diff(_xi[r](t, r, S, V), t))*r*S+4*_xi[S](t, r, S, V)*`λΣ`*r^(1+beta)-4*r*(a-`κr`)*_xi[S](t, r, S, V)-2*beta*`λΣ`*r^beta*_xi[r](t, r, S, V)*S, diff(diff(_xi[S](t, r, S, V), V), V), diff(diff(_xi[r](t, r, S, V), V), V), diff(diff(_xi[t](t, r, S, V), V), V), diff(diff(_xi[t](t, r, S, V), V), r), diff(_xi[S](t, r, S, V), V), diff(_xi[r](t, r, S, V), V), diff(_xi[t](t, r, S, V), S), diff(_xi[t](t, r, S, V), V), diff(_xi[t](t, r, S, V), r)}

(3)

Try now processing this taking beta as an independent variable, that simplifies the problem significantly

PDEtools:-casesplit({(S*(diff(_xi[r](t, r, S, V), r))-S*(diff(_xi[S](t, r, S, V), S))+_xi[S](t, r, S, V))*r^(1+2*beta)-beta*r^(2*beta)*_xi[r](t, r, S, V)*S, (diff(diff(_xi[t](t, r, S, V), S), V))*S^2*sigma^2-2*(diff(_xi[S](t, r, S, V), V)), -2*r^(1+2*beta)*(diff(_xi[S](t, r, S, V), r))*S*Sigma^2-2*(diff(_xi[r](t, r, S, V), S))*r*S^3*sigma^2, diff(diff(_xi[r](t, r, S, V), V), r)-(1/2)*(diff(diff(_eta[V](t, r, S, V), V), V)), r^(1+2*beta)*(diff(diff(_xi[S](t, r, S, V), V), r))*Sigma^2+(diff(diff(_xi[r](t, r, S, V), S), V))*S^2*r*sigma^2+2*(`λΣ`*r^(1+beta)-r*(a-`κr`))*(diff(_xi[S](t, r, S, V), V)), (diff(diff(_xi[S](t, r, S, V), S), V))*S*sigma^2-(1/2)*(diff(diff(_eta[V](t, r, S, V), V), V))*S*sigma^2-2*(diff(_xi[S](t, r, S, V), V)), r^(1+2*beta)*(diff(diff(_eta[V](t, r, S, V), r), r))*S*Sigma^2+S^3*r*(diff(diff(_eta[V](t, r, S, V), S), S))*sigma^2-2*S*(`λΣ`*r^(1+beta)-r*(a-`κr`))*(diff(_eta[V](t, r, S, V), r))+2*(S*r*V*(diff(_eta[V](t, r, S, V), V))+S^2*(diff(_eta[V](t, r, S, V), S))+S*(diff(_eta[V](t, r, S, V), t))-2*S*(diff(_xi[S](t, r, S, V), S))*r*V-S*r*_eta[V](t, r, S, V)-V*(S*_xi[r](t, r, S, V)-2*_xi[S](t, r, S, V)*r))*r, -(diff(diff(_xi[t](t, r, S, V), r), r))*r^(1+2*beta)*S*Sigma^2-S^3*r*(diff(diff(_xi[t](t, r, S, V), S), S))*sigma^2+2*S*(`λΣ`*r^(1+beta)-r*(a-`κr`))*(diff(_xi[t](t, r, S, V), r))-2*r*(S*(diff(_xi[t](t, r, S, V), V))*V*r+S^2*(diff(_xi[t](t, r, S, V), S))+S*(diff(_xi[t](t, r, S, V), t))-2*S*(diff(_xi[S](t, r, S, V), S))+2*_xi[S](t, r, S, V)), r^(1+2*beta)*(diff(diff(_xi[S](t, r, S, V), r), r))*Sigma^2-2*(diff(diff(_eta[V](t, r, S, V), S), V))*S^2*r*sigma^2+(diff(diff(_xi[S](t, r, S, V), S), S))*S^2*r*sigma^2+(-2*`λΣ`*r^(1+beta)+2*r*(a-`κr`))*(diff(_xi[S](t, r, S, V), r))-2*r*(-3*(diff(_xi[S](t, r, S, V), V))*V*r+S*(diff(_xi[S](t, r, S, V), S))-_xi[S](t, r, S, V)-(diff(_xi[S](t, r, S, V), t))), -S*Sigma^2*(diff(diff(_xi[r](t, r, S, V), r), r)-2*(diff(diff(_eta[V](t, r, S, V), V), r)))*r^(1+2*beta)-r*(diff(diff(_xi[r](t, r, S, V), S), S))*S^3*sigma^2+2*(`λΣ`*r^(1+beta)-r*(a-`κr`))*S*(diff(_xi[r](t, r, S, V), r))-4*(`λΣ`*r^(1+beta)-r*(a-`κr`))*S*(diff(_xi[S](t, r, S, V), S))-2*(diff(_xi[r](t, r, S, V), V))*V*r^2*S-2*(diff(_xi[r](t, r, S, V), S))*r*S^2-2*(diff(_xi[r](t, r, S, V), t))*r*S+4*_xi[S](t, r, S, V)*`λΣ`*r^(1+beta)-4*r*(a-`κr`)*_xi[S](t, r, S, V)-2*beta*`λΣ`*r^beta*_xi[r](t, r, S, V)*S, diff(diff(_xi[S](t, r, S, V), V), V), diff(diff(_xi[r](t, r, S, V), V), V), diff(diff(_xi[t](t, r, S, V), V), V), diff(diff(_xi[t](t, r, S, V), V), r), diff(_xi[S](t, r, S, V), V), diff(_xi[r](t, r, S, V), V), diff(_xi[t](t, r, S, V), S), diff(_xi[t](t, r, S, V), V), diff(_xi[t](t, r, S, V), r)}, ivars = [t, r, S, V, beta])

`casesplit/ans`([diff(_eta[V](t, r, S, V), t) = (1/2)*((diff(_xi[S](t, r, S, V), t))*V*sigma^2-2*(diff(_xi[S](t, r, S, V), t))*V)/(sigma^2*S), diff(_eta[V](t, r, S, V), r) = 0, diff(_eta[V](t, r, S, V), S) = (diff(_xi[S](t, r, S, V), t))*V/(sigma^2*S^2), diff(_eta[V](t, r, S, V), V) = _eta[V](t, r, S, V)/V, diff(diff(_xi[S](t, r, S, V), t), t) = 0, diff(_xi[S](t, r, S, V), r) = 0, diff(_xi[S](t, r, S, V), S) = _xi[S](t, r, S, V)/S, diff(_xi[S](t, r, S, V), V) = 0, _xi[r](t, r, S, V) = 0, diff(_xi[t](t, r, S, V), t) = 0, diff(_xi[t](t, r, S, V), r) = 0, diff(_xi[t](t, r, S, V), S) = 0, diff(_xi[t](t, r, S, V), V) = 0], [])

(4)

Solve it now to compute the infinitesimals symmety generators themselves

pdsolve(`casesplit/ans`([diff(_eta[V](t, r, S, V), t) = (1/2)*((diff(_xi[S](t, r, S, V), t))*V*sigma^2-2*(diff(_xi[S](t, r, S, V), t))*V)/(sigma^2*S), diff(_eta[V](t, r, S, V), r) = 0, diff(_eta[V](t, r, S, V), S) = (diff(_xi[S](t, r, S, V), t))*V/(sigma^2*S^2), diff(_eta[V](t, r, S, V), V) = _eta[V](t, r, S, V)/V, diff(diff(_xi[S](t, r, S, V), t), t) = 0, diff(_xi[S](t, r, S, V), r) = 0, diff(_xi[S](t, r, S, V), S) = _xi[S](t, r, S, V)/S, diff(_xi[S](t, r, S, V), V) = 0, _xi[r](t, r, S, V) = 0, diff(_xi[t](t, r, S, V), t) = 0, diff(_xi[t](t, r, S, V), r) = 0, diff(_xi[t](t, r, S, V), S) = 0, diff(_xi[t](t, r, S, V), V) = 0], []))

{_eta[V](t, r, S, V) = (1/2)*V*(2*_C1*ln(S)+(_C1*t+2*_C4)*sigma^2-2*_C1*t)/sigma^2, _xi[S](t, r, S, V) = (_C1*t+_C2)*S, _xi[r](t, r, S, V) = 0, _xi[t](t, r, S, V) = _C3}

(5)

Verify that this indeed solves the Determining system for the symmetries= (note that I test directly against (3), not the simplified form (4), that is also canceled by (5)):

pdetest({_eta[V](t, r, S, V) = (1/2)*V*(2*_C1*ln(S)+(_C1*t+2*_C4)*sigma^2-2*_C1*t)/sigma^2, _xi[S](t, r, S, V) = (_C1*t+_C2)*S, _xi[r](t, r, S, V) = 0, _xi[t](t, r, S, V) = _C3}, {(S*(diff(_xi[r](t, r, S, V), r))-S*(diff(_xi[S](t, r, S, V), S))+_xi[S](t, r, S, V))*r^(1+2*beta)-beta*r^(2*beta)*_xi[r](t, r, S, V)*S, (diff(diff(_xi[t](t, r, S, V), S), V))*S^2*sigma^2-2*(diff(_xi[S](t, r, S, V), V)), -2*r^(1+2*beta)*(diff(_xi[S](t, r, S, V), r))*S*Sigma^2-2*(diff(_xi[r](t, r, S, V), S))*r*S^3*sigma^2, diff(diff(_xi[r](t, r, S, V), V), r)-(1/2)*(diff(diff(_eta[V](t, r, S, V), V), V)), r^(1+2*beta)*(diff(diff(_xi[S](t, r, S, V), V), r))*Sigma^2+(diff(diff(_xi[r](t, r, S, V), S), V))*S^2*r*sigma^2+2*(`λΣ`*r^(1+beta)-r*(a-`κr`))*(diff(_xi[S](t, r, S, V), V)), (diff(diff(_xi[S](t, r, S, V), S), V))*S*sigma^2-(1/2)*(diff(diff(_eta[V](t, r, S, V), V), V))*S*sigma^2-2*(diff(_xi[S](t, r, S, V), V)), r^(1+2*beta)*(diff(diff(_eta[V](t, r, S, V), r), r))*S*Sigma^2+S^3*r*(diff(diff(_eta[V](t, r, S, V), S), S))*sigma^2-2*S*(`λΣ`*r^(1+beta)-r*(a-`κr`))*(diff(_eta[V](t, r, S, V), r))+2*(S*r*V*(diff(_eta[V](t, r, S, V), V))+S^2*(diff(_eta[V](t, r, S, V), S))+S*(diff(_eta[V](t, r, S, V), t))-2*S*(diff(_xi[S](t, r, S, V), S))*r*V-S*r*_eta[V](t, r, S, V)-V*(S*_xi[r](t, r, S, V)-2*_xi[S](t, r, S, V)*r))*r, -(diff(diff(_xi[t](t, r, S, V), r), r))*r^(1+2*beta)*S*Sigma^2-S^3*r*(diff(diff(_xi[t](t, r, S, V), S), S))*sigma^2+2*S*(`λΣ`*r^(1+beta)-r*(a-`κr`))*(diff(_xi[t](t, r, S, V), r))-2*r*(S*(diff(_xi[t](t, r, S, V), V))*V*r+S^2*(diff(_xi[t](t, r, S, V), S))+S*(diff(_xi[t](t, r, S, V), t))-2*S*(diff(_xi[S](t, r, S, V), S))+2*_xi[S](t, r, S, V)), r^(1+2*beta)*(diff(diff(_xi[S](t, r, S, V), r), r))*Sigma^2-2*(diff(diff(_eta[V](t, r, S, V), S), V))*S^2*r*sigma^2+(diff(diff(_xi[S](t, r, S, V), S), S))*S^2*r*sigma^2+(-2*`λΣ`*r^(1+beta)+2*r*(a-`κr`))*(diff(_xi[S](t, r, S, V), r))-2*r*(-3*(diff(_xi[S](t, r, S, V), V))*V*r+S*(diff(_xi[S](t, r, S, V), S))-_xi[S](t, r, S, V)-(diff(_xi[S](t, r, S, V), t))), -S*Sigma^2*(diff(diff(_xi[r](t, r, S, V), r), r)-2*(diff(diff(_eta[V](t, r, S, V), V), r)))*r^(1+2*beta)-r*(diff(diff(_xi[r](t, r, S, V), S), S))*S^3*sigma^2+2*(`λΣ`*r^(1+beta)-r*(a-`κr`))*S*(diff(_xi[r](t, r, S, V), r))-4*(`λΣ`*r^(1+beta)-r*(a-`κr`))*S*(diff(_xi[S](t, r, S, V), S))-2*(diff(_xi[r](t, r, S, V), V))*V*r^2*S-2*(diff(_xi[r](t, r, S, V), S))*r*S^2-2*(diff(_xi[r](t, r, S, V), t))*r*S+4*_xi[S](t, r, S, V)*`λΣ`*r^(1+beta)-4*r*(a-`κr`)*_xi[S](t, r, S, V)-2*beta*`λΣ`*r^beta*_xi[r](t, r, S, V)*S, diff(diff(_xi[S](t, r, S, V), V), V), diff(diff(_xi[r](t, r, S, V), V), V), diff(diff(_xi[t](t, r, S, V), V), V), diff(diff(_xi[t](t, r, S, V), V), r), diff(_xi[S](t, r, S, V), V), diff(_xi[r](t, r, S, V), V), diff(_xi[t](t, r, S, V), S), diff(_xi[t](t, r, S, V), V), diff(_xi[t](t, r, S, V), r)})

{0}

(6)

So, for an infinitesimal of this form,

[_xi[t](t, r, S, V), _xi[r](t, r, S, V), _xi[S](t, r, S, V), _eta[V](t, r, S, V)]

[_xi[t](t, r, S, V), _xi[r](t, r, S, V), _xi[S](t, r, S, V), _eta[V](t, r, S, V)]

(7)

This is a symmetry group involving three arbitrary constants

X := eval([_xi[t](t, r, S, V), _xi[r](t, r, S, V), _xi[S](t, r, S, V), _eta[V](t, r, S, V)], {_eta[V](t, r, S, V) = (1/2)*V*(2*_C1*ln(S)+(_C1*t+2*_C4)*sigma^2-2*_C1*t)/sigma^2, _xi[S](t, r, S, V) = (_C1*t+_C2)*S, _xi[r](t, r, S, V) = 0, _xi[t](t, r, S, V) = _C3})

[_C3, 0, (_C1*t+_C2)*S, (1/2)*V*(2*_C1*ln(S)+(_C1*t+2*_C4)*sigma^2-2*_C1*t)/sigma^2]

(8)

Alternatively you can also verify that this is indeed a symmetry directly against the PDE equ

PDEtools:-SymmetryTest(X, equ)

{0}

(9)

Time to compute solutions. You can do that without specializing the symmetry, arriving at a solution involving Airy functions and the solutions of a 2nd order linear ODE that the system fails to solve (you'd need to specialize beta to obtain something more concrete)

PDEtools:-InvariantSolutions(equ, HINT = X)

V(t, r, S) = DESol({(-2*r^(-2*beta+1)*_Y(r)+((2*a-2*`κr`)*(diff(_Y(r), r))-_Y(r)*_c[1])*r^(-2*beta)-2*r^(-beta)*(diff(_Y(r), r))*`λΣ`+(diff(diff(_Y(r), r), r))*Sigma^2)/Sigma^2}, {_Y(r)})*(S*exp(-(1/2)*t*(_C1*t+2*_C2)/_C3))^((1/2)*(_C8*sigma^2+2*_C7-2*_C8)/(sigma^2*_C8))*(AiryAi((1/8)*2^(1/3)*(-8*ln(S*exp(-(1/2)*t*(_C1*t+2*_C2)/_C3))*_C6*_C8+(4+sigma^4+(-4*_c[1]-4)*sigma^2)*_C8^2+((4*_C7-8*_C9)*sigma^2-8*_C7)*_C8+4*_C7^2)*(_C6/(_C8*sigma^4))^(1/3)/(_C6*_C8))*_C10+AiryBi((1/8)*2^(1/3)*(-8*ln(S*exp(-(1/2)*t*(_C1*t+2*_C2)/_C3))*_C6*_C8+(4+sigma^4+(-4*_c[1]-4)*sigma^2)*_C8^2+((4*_C7-8*_C9)*sigma^2-8*_C7)*_C8+4*_C7^2)*(_C6/(_C8*sigma^4))^(1/3)/(_C6*_C8))*_C5)/(exp((1/3)*t*(_C1^2*t^2+(3/2)*_C1*(-(1/2)*sigma^2*_C3+_C2+_C3)*t-3*_C3*_C4*sigma^2)/(_C3^2*sigma^2))*S^(-_C1*t/(_C3*sigma^2)))

(10)

Typically (not in this case), specializing the symmetry helps, and for this purpose you can use the PDEtools:-Library:- Specialize_Cn command  (see ? PDEtools:-Library )

Y := PDEtools:-Library:-Specialize_Cn(X)

[0, 0, t*S, (1/2)*V*(2*ln(S)+t*sigma^2-2*t)/sigma^2], [0, 0, S, 0], [1, 0, 0, 0], [0, 0, 0, V]

(11)

Verify these symmetries

map(PDEtools:-SymmetryTest, [[0, 0, t*S, (1/2)*V*(2*ln(S)+t*sigma^2-2*t)/sigma^2], [0, 0, S, 0], [1, 0, 0, 0], [0, 0, 0, V]], equ)

[{0}, {0}, {0}, {0}]

(12)

You can now use each of these to generate different kinds of invariant solutions, although as said in this example all of them lead to solutions to equ that depend on the solution of a 2nd order linear ODE that the system fails to solve (you'd need to specialize beta to arrive somewhere more concrete). For example:

PDEtools:-InvariantSolutions(equ, HINT = Y[1])

V(t, r, S) = _C1*exp((1/2)*_c[1]*t)*DESol({(1/4)*(-8*r^(-2*beta+1)*_Y(r)*sigma^2+(8*sigma^2*(a-`κr`)*(diff(_Y(r), r))-(4+sigma^4+(-4*_c[1]-4)*sigma^2)*_Y(r))*r^(-2*beta)-8*(diff(_Y(r), r))*r^(-beta)*sigma^2*`λΣ`+4*(diff(diff(_Y(r), r), r))*Sigma^2*sigma^2)/(Sigma^2*sigma^2)}, {_Y(r)})/(t^(1/2)*S^(-(1/2)*(sigma^2-2)/sigma^2)*exp(-(1/2)*ln(S)^2/(t*sigma^2)))

(13)

PDEtools:-InvariantSolutions(equ, HINT = Y[2])

V(t, r, S) = _C1*exp(_c[1]*t)*DESol({(-2*r^(-2*beta+1)*_Y(r)+((2*a-2*`κr`)*(diff(_Y(r), r))+2*_Y(r)*_c[1])*r^(-2*beta)-2*r^(-beta)*(diff(_Y(r), r))*`λΣ`+(diff(diff(_Y(r), r), r))*Sigma^2)/Sigma^2}, {_Y(r)})

(14)

PDEtools:-InvariantSolutions(equ, HINT = Y[3])

V(t, r, S) = S^(1/2)*DESol({(-2*r^(-2*beta+1)*_Y(r)+((2*a-2*`κr`)*(diff(_Y(r), r))-_Y(r)*_c[1])*r^(-2*beta)-2*r^(-beta)*(diff(_Y(r), r))*`λΣ`+(diff(diff(_Y(r), r), r))*Sigma^2)/Sigma^2}, {_Y(r)})*(S^((1/2)*(-2+(4+sigma^4+(-4*_c[1]-4)*sigma^2)^(1/2))/sigma^2)*_C1+S^(-(1/2)*(2+(4+sigma^4+(-4*_c[1]-4)*sigma^2)^(1/2))/sigma^2)*_C2)

(15)

"And no invariant solution is at reach for Y[4]."

 

You can also play around directly with the options of InvariantSolutions. For example,

PDEtools:-InvariantSolutions(equ, degreeofinfinitesimals = 1)

V(t, r, S) = DESol({(-2*r^(-2*beta+1)*_Y(r)+2*(diff(_Y(r), r))*(a-`κr`)*r^(-2*beta)-2*r^(-beta)*(diff(_Y(r), r))*`λΣ`+(diff(diff(_Y(r), r), r))*Sigma^2)/Sigma^2}, {_Y(r)})*_C1

(16)

pdetest(%, equ)

0

(17)

But in this example I think it is improbable that you will obtain an actual solution to equ without specializing the constants.


Download Equ_InvariantSolutions.mw

Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions, Maplesoft

Hi
The conversion network for mathematical functions in Maple is wonderful/powerful. You can take advantage of it for this problem. Besides Carl's solution, you can also try directly a conversion to elementary form, for instance via

> conert(...<you hypergeometric expression here> ... , elementary)

Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions, Maplesoft

The PDEtools package has a complete set of commands for computing symmetries of partial differential equations or systems of them (31 commands). The determinant system for the symmetries are computed using PDEtools:-DeterminantPDE, then you can solve that system using pdsolve. But, simpler, you can compute the symmetries directly, using PDEtools:-Infinitesimals. All this is explained in the help page ?PDEtools.

Edgardo S. Cheb-Terrab
Physics, Maplesoft

 

Hi Glowing

 

Interesting example yours. First: given an ODE & ICs, the "mandate" of dsolve is to find one solution that solves the problem, assuming, generally speaking, that this solution is unique, because, generally speaking, an nth order ODE with n initial/boundary conditions has only 1 solution.

 

There are examples, however, where the solution is not unique, there is more than one solution, the additional ones typically associated with the existence of singularities of the ODE, or branches (e.g. sqrt in your example).

 

What you can expect from dsolve in these more tricky cases, where the system of an nth order ODE & n ICS has more than one solution, is one solution only, the one dsolve can find faster. This is so because dsolve counts the initial conditions, checks the differential order, and if these numbers are the same it automatically runs for the simpler solution to the problem, as said, assuming that in these cases the solution is unique...

 

In this framework, the solution returned U(t) = 0 is correct, it cancels the ODE and the ics. Nothing wrong with that.

 

But then: how do you compute the other solutions? Although in general you do not need to use optional arguments, in cases like this one you need: use the 'usesolutions = ...'. This optional argument is explained in the help page ?dsolve,details, the section on Optional Arguments. The possible right-hand sides are explained in the help page and to get the solution you want you can use "general", as in

dsolve({ics, ode}, usesolutions = "general");

U(t) = U__0*((exp(t))^2-2*exp(t)+1)/(exp(t))^2

(1)

This is the same solution you derived step-by-step. To get both solutions use "particular via symmetries and general", as in

dsolve({ics, ode}, usesolutions = "particular via symmetries and general");

U(t) = 0, U(t) = U__0*(exp(t)-1)^2/(exp(t))^2

(2)

So up to here: the solution returned by dsolve without optional arguments is correct, it cancels the ODE and ics, so it is not a hidden bug; the problem passed is special, admits more than one solution, to get the other solution you can use dsolve's optional argument 'usesolutions'.

 

There is still one things: how can you move away from this situation where the set of a nth order ODE & n ICs has more than the one unique solution? Generally speaking, the answer is to avoid constructing the ODE through differentiation of expressions with branches, as for instance sqrt(U(t)) is. In your example you constructed your ODE via

ode := diff(sqrt(U(t)), t) = sqrt(U__0)-sqrt(U(t))

You can verify the existence of branches using FunctionAdvisor(branch_cuts, expression), say as in

FunctionAdvisor(branch_cuts, sqrt(U))

[U^(1/2), U < 0]

(3)

You see that the cut starts at 0, precisely where you gave your initial condition. Not a good a idea. So, in this example you posted an easy manner to move away from this situation is to replace, before constructing the ODE through differentiation, any occurrence of sqrt(U(t)) by - say - V(t), then U(0) = 0 becomes V(0) = 0, then the system ODE + ICS becomes well defined and has a unique solution - corresponding to the one you were expecting. In your example, constructing the initial condition for V is trivial; in non-trivial situations you can use dchange, for instance for your example it gives

PDEtools:-dchange({t = T, U(t) = V(T)}, U(0))

V(0)

(4)

In summary: all OK with the solution by dsolve without using optional arguments and you can use the usesolutions option to get the other solution. You can avoid these situations by not constructing the ODE problem through differentiation of an expression with a branch cut while giving the initial condition at the origin or over the cut.

 

Edgardo S. Cheb-Terrab

Physics, Maplesoft

 

 

Download Option_usesolution.mw

 

Hi Mac
In brief, right after 'with(Physics:-Vectors)' enter 'Setup(geometricdifferentiation = false)', and then the flow will be as you expected (already verified that with the worksheet you posted). For details see ?Physics:-Setup, where geometricdifferentiation is explained, or with more details see ?Physics[Vectors][diff].

Alternatively, you can perform the same computation without canceling geometricdifferentiation: you'd need to indicate to dchange also the change rule for the dependent variable of your problem (convenient in general), that is E1_(x,y,z,t), so your input line for dchange would be:

TR := {x = r*cos(theta), z = r*sin(theta), E1_(x,y,z,t) = R(r)*exp(I*(-k*y+omega*t))};

PDEtools:-dchange(TR, -(diff(E1_(x, y, z, t), x, x))-(diff(E1_(x, y, z, t), y, y))-(diff(E1_(x, y, z, t), z, z)) = -mu*epsilon*(diff(E1_(x, y, z, t), t, t)), [r, theta, R(r)], simplify);

This is a better approach, keeps geometridifferentiation working, and removes the need of that posterior subs that you do in the worksheet.

Regarding updating Physics in Maple 17.02, have no doubt, it is convenient: you will have access to many improvements as well as most of Maple 18 developments, and bug fixes. The place were you get the update is always the same: the Maplesoft R&D page for Physics. There you see a link to the Maple 18 updates that happened so far (the library in this link does not work for Maple 17) and another link for the latest update for Maple 17, dated March/18/2014.

Edgardo S. Cheb-Terrab
Physics, Maplesoft

Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions, Maplesoft

Hi Mac,
This is yet another solution for Maples before Maple 18:

> unprotect(intat):
> Intat := setattribute('Intat', '_syslib'):
> protect(intat):

That will just make Intat fully inert and resolve your issue the same way as the other 'workarounds', though in this case it is more a solution than a workaround. And for Maple 18, this will be available in a couple of days on the webpage with updates to the DE and mathematical functions libraries.

Edgardo S. Cheb-Terrab
Physics, Maplesoft

Hi
Indeed, Alejandro is right when saying that Intat is not entirely inert: it checks arguments, from a time (actually, 1996) when inert routines were under discussion. I will remove this check of arguments.

Then you can also use %intat (so: the active lowercase command preceeded by the %) instead of Intat, and if that doesn't work, it needs to be fixed. I have not tried with your example, but am cooking an update on DEs and Mathematical Functions to be released this week and can include a fix for this too if it doesn't work (please let me know).

Then there is the use of single versus double quoting or use of value: generally speaking, double quoting or double calls to value should not be necessary, but if you use %Intat (so the uppercase preceeded by %) you will naturally need double calls to value, the same way if you use %Int instead of %int, for example as in value(value(%Int(x,x))).

The message I am trying to pass is that the real issues are to be fixed, not worked around. The fixes are easy to accomplish and can be available to everybody as soon as this week (CAVEAT: they only work in Maple 18).

Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions, Maplesoft

Hi
Just gave a look and updated the Maplesoft webpage: this issue related to having to apply Simplify twice to obtain zero is fixed. So now Simplify((16) - (14)) returns zero in one go and the same for Simplify((rhs((17)) - rhs((15))).

Edgardo S. Cheb-Terrab
Physics, Maplesoft

Hi Markiyan

In cases like this one (fractional powers and some other ones) one approach is to transform your problem into one that is a bit more general, solve that one, then try restricting the solutions to solve the original problem - sometimes you only need to select the solution that is useful.

For example, you starting expression is

y^2+y^3+(y^3-x^2-3*x*y)^(1/4) <= 5*x*y;

Isolate the radical

isolate(%, (y^3-x^2-3*x*y)^(1/4));

Now map both sides to the power 4

map(u -> u^4, %);

Solve this more-general problem:

solve(%, allsolutions);  # bunch of solutions here

And, I haven't tried, but either one of these already solves your problem or from their remainder (obtained after substituting the solution into the original problem and simplifying) you could see how to particularize it further.

Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions, Maplesoft

Looking closer, the simplification works fine. You say that "equation (16) should correspond to equation (14)" and that is indeed the case. One way to verify that is to take the difference (16) (14) and perform the summation over all the repeated indices:

Vacuum_Solutions_reviewed.mw

Edgardo S. Cheb-Terrab
Physics, Maplesoft

Hi Mac,
There is no addBasis command, the three coordinate systems with which the Physics :-Vectors package works are fixed, as explained in ?Physics,Vectors : [x, y, z], [rho, phi, z], [r, theta, phi]. I do have in mind implementing arbitrary basis, but it is not high in the list of priorities, mainly because these three systems of coordinates cover most applications & textbooks, swapping two coordinates in the basis definition does not really change the geometry (the existing system is mostly as good for that purpose), and there are various other things - new functionality - that look more relevant or having higher positive impact if implemented.

Regarding Physics being open: it is as open as any other Maple program. You can see the routines and change them if you want. More important: Physics (starting with Maple 17) comes with a library of more than 100 dedicated-physics-programming commands, plus approx 70 physics-dedicated-types, basically all the routines used to write Physics, that allow you to extend the package in different ways if you are willing to program. Various of these "programming" routines are also useful at interactive level. Details are in the help page (available online) ?Physics,Library

But you need to update your Maple, Mac. Have in mind that the Physics package duplicated its size for Maple 16, then again for Maple 17, and mostly one more time for Maple 18, not to mention that by now more than 40 updates were posted in the Maplesoft R&D Physics webpage after Maple 18 got released.

Edgardo S. Cheb-Terrab
Physics, Maplesoft

Hi Mac,

The problem in old maples is that assuming replaces the cartesian coordiante x with a local one and then the Nabla operator in Physics differentiates as in diff(local_x, global_x) = 0. This has been fixed in previous Maples - it works fine in Maple 18, and also in Maple 17 + the Physics update. This is what I receive in Maple 18:

 

"Maple Initialization loaded..."

(1)

restart:with(Physics[Vectors]);
Setup(mathematicalnotation = true);

[`&x`, `+`, `.`, ChangeBasis, Component, Curl, DirectionalDiff, Divergence, Gradient, Identify, Laplacian, Nabla, Norm, Setup, diff]

 

[mathematicalnotation = true]

(2)

# Maxwell's eqn
M4 := `&x`(Nabla, B1_(x, y, z, t)) = mu*epsilon*(diff(E1_(x, y, z, t), t));

Physics:-Vectors:-Curl(B1_(x, y, z, t)) = mu*epsilon*(diff(E1_(x, y, z, t), t))

(3)

eval(subs(B1_(x,y,z,t)=Bxx(x,y,z,t)*_i+Bzz(x,y,z,t)*_k,M4));

(diff(Bzz(x, y, z, t), y))*_i+(diff(Bxx(x, y, z, t), z)-(diff(Bzz(x, y, z, t), x)))*_j-(diff(Bxx(x, y, z, t), y))*_k = mu*epsilon*(diff(E1_(x, y, z, t), t))

(4)

eval(subs(B1_(x,y,z,t)=Bxx(x,y,z,t)*_i+Bzz(x,y,z,t)*_k,M4)) assuming real;

(diff(Bzz(x, y, z, t), y))*_i+(diff(Bxx(x, y, z, t), z)-(diff(Bzz(x, y, z, t), x)))*_j-(diff(Bxx(x, y, z, t), y))*_k = mu*epsilon*(diff(E1_(x, y, z, t), t))

(5)

The above is in Maple 18, but you can reproduce it in Maple 17 as well including the update of Physics available on the Maplesoft R&D Physics webpage.

 

In brief: you need to update your Maple.

 

NULL

 

Download Maxwell_testMaple18.mw

Edgardo S. Cheb-Terrab
Physics, Maplesoft

Hi

I am not sure I understand your problem, but if what you want is to solve 

ee := A*cos(x)+R*i*sin(x) = B*exp(-x)+C*exp(x);

For A, B, C, R independent of x, and such that ee is true for arbitrary values of x, as far as I can tell the problem has no solution but for the trivial solution. The command in Maple that performs this computation is PDEtools:-Solve. Check its help page for details. The input/output for your example is as follows:

PDEtools:-Solve(ee, {A, B, C, R}, independentof = x);

                                                  {A = 0, B = 0, C = 0, R = 0}


Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions, Maplesoft

@John Fredsted 

We included the fix to this problem in the latest update for Maple 17, that also includes all the changes in the package till March 18 2014, and with this we close the updates for Maple 17. This update is available for download as usuai on the Maplesoft R&D Physics webpage (note that there are two links: one for the ongoing updates in Maple 18, these do not work in Maple 17, and another link, Physics-61.3.mla.zip, with the latest update for Maple 17).

Edgardo S. Cheb-Terrab
Physics, Differential Equations and Mathematical Functions, Maplesoft.

First 48 49 50 51 52 53 54 Last Page 50 of 60