Advanced Search

nm

9934 Reputation

19 Badges

12 years, 90 days
Contact nm

MaplePrimes Activity


These are answers submitted by nm

reply...

nm 9934
October 09 2024
1 2

set the ranges then it works like Matlab

> 

f:=x^(2*z)+y^(-2*y^(-z))+exp(-1/10*z^2)=1;
plots:-implicitplot3d(f,x=1.3..2,y=1.3..2,z=-2..-1.3,style=surface)
> 

 

x^(2*z)+y^(-2*y^(-z))+exp(-(1/10)*z^2) = 1

> 

 


As for the gridlines, this option does not seem to be supported for 3D plots in Maple. May be someone knows of a trick.

 

Download plot3d.mw

reply...

nm 9934
October 08 2024
0 0

Root of is just symbolic representation of roots of equations. Many cas system use it. 

To remove it, try calling allvalues() on the result. In your case, it gives 

             `[Length of output exceeds limit of 1000000]`

So it is very large result and something you probably do not want to see, so better keep the RoofOf there  :)

reply...

nm 9934
October 06 2024
1 1

This fixes two issues you had. you need to remove the O() from the series solution, and need to add [] around ode and the IC. But Maple can't solve the final pde analytically.  V 2024.1
 

> 

restart;
> 

ode0 := diff(xi^2*diff(theta[E](xi), xi), xi)/xi^2 = -theta[E](xi)^n;

(2*xi*(diff(theta[E](xi), xi))+xi^2*(diff(diff(theta[E](xi), xi), xi)))/xi^2 = -theta[E](xi)^n

> 

bc0 := theta[E](0) = 1, D(theta[E])(0) = 0;

theta[E](0) = 1, (D(theta[E]))(0) = 0

> 

base := dsolve([bc0, ode0], theta[E](xi), series);
base := convert(base,polynom);

theta[E](xi) = series(1-(1/6)*xi^2+((1/120)*n)*xi^4+O(xi^6),xi,6)

theta[E](xi) = 1-(1/6)*xi^2+(1/120)*n*xi^4

> 

pde1 := diff(xi^2*diff(psi(xi, mu), xi), xi)/xi^2 + diff((-mu^2 + 1)*diff(psi(xi, mu), mu), mu)/xi^2 = -psi(xi, mu)^n + v;

(2*xi*(diff(psi(xi, mu), xi))+xi^2*(diff(diff(psi(xi, mu), xi), xi)))/xi^2+(-2*mu*(diff(psi(xi, mu), mu))+(-mu^2+1)*(diff(diff(psi(xi, mu), mu), mu)))/xi^2 = -psi(xi, mu)^n+v

> 

bc1 := psi(0, mu) = 1, D[1](psi)(0, mu) = 0, D[2](psi)(0, mu) = 0, limit(psi(xi, mu), v = 0) = rhs(base);

psi(0, mu) = 1, (D[1](psi))(0, mu) = 0, (D[2](psi))(0, mu) = 0, psi(xi, mu) = 1-(1/6)*xi^2+(1/120)*n*xi^4

> 

psi(xi, mu)

psi(xi, mu)

> 

pdsolve([pde1, bc1],psi(xi, mu))
> 

pdsolve([pde1, bc1],psi(xi, mu),series)
> 

pdsolve(pde1,psi(xi, mu))
> 

 


 

Download Nonlinear_Elliptic_PDE_in_Spherical_Coordinate.mw

reply...

nm 9934
October 05 2024
1 0

Not good in using paint.exe to draw triangles. This is the best I could do. In this diagram, A is the length of one side of the largest square inside the big triangle. (there can only be one such square).

Area of big triangle = area of (1) + area of (2) + area of (3) + area of inside square.

But area of big triangle = 1/2*h*c

area of (1) = 1/2*A*c1 

area of (2) = 1/2*(h-A)*A

area of (3) = 1/2*A*c2

And we have that c1+A+c2 = C

Set up equation and solve for A, this gives A= (h*c)/(c+h).

Now we do not know h, the height of the main triangle. To get rid of h, we use heron's formula which finds area of triangle using perimeter of triangle. Here is updated worksheet. The answer gives A as function of the unknown C side of the triangle.

Then use Optimization:-Maximize to maximize the above for c>=4 (since you said the two smaller sides are 3 and 4) and Maple gave A=1.65425211895987911 


 

> 

restart;
> 

eq1:=1/2*h*c= 1/2*A*c1 + 1/2*A*c2 + 1/2*(h-A)*A + A^2;
eq2:=c1+c2+A=c;

(1/2)*h*c = (1/2)*A*c1+(1/2)*A*c2+(1/2)*(h-A)*A+A^2

c1+c2+A = c

> 

sol:=solve([eq1,eq2],[A,c],explicit)

[[A = -(1/2)*c1-(1/2)*c2+(1/2)*(c1^2+2*c1*c2+4*c1*h+c2^2+4*c2*h)^(1/2), c = (1/2)*c1+(1/2)*c2+(1/2)*(c1^2+2*c1*c2+4*c1*h+c2^2+4*c2*h)^(1/2)], [A = -(1/2)*c1-(1/2)*c2-(1/2)*(c1^2+2*c1*c2+4*c1*h+c2^2+4*c2*h)^(1/2), c = (1/2)*c1+(1/2)*c2-(1/2)*(c1^2+2*c1*c2+4*c1*h+c2^2+4*c2*h)^(1/2)]]

> 

eq:=lhs(sol[1,1])=simplify(rhs(sol[1,1]),[c1+c2+A=c])

A = (1/2)*((A-c)*(A-c-4*h))^(1/2)-(1/2)*c+(1/2)*A

> 

side_of_square:=PDEtools:-Solve(eq,A)

A = h*c/(c+h)

> 

#Use heron's formula to find h
a:=3;b:=4;
s:=(a+b+c)/2; #semi perimeter
area_of_big_triangle:= sqrt(s*(s-a)*(s-b)*(s-c));

3

4

7/2+(1/2)*c

((7/2+(1/2)*c)*(1/2+(1/2)*c)*(-1/2+(1/2)*c)*(7/2-(1/2)*c))^(1/2)

> 

eq2:=1/2*h*c = area_of_big_triangle

(1/2)*h*c = ((7/2+(1/2)*c)*(1/2+(1/2)*c)*(-1/2+(1/2)*c)*(7/2-(1/2)*c))^(1/2)

> 

h_sol:=PDEtools:-Solve(eq2,h);

h = (1/2)*(-(c^2-1)*(c^2-49))^(1/2)/c

> 

side_of_square:=simplify(eval(side_of_square,h_sol));

A = (-c^4+50*c^2-49)^(1/2)*c/(2*c^2+(-c^4+50*c^2-49)^(1/2))

> 

side_of_square:=Optimization:-Maximize(rhs(side_of_square),{c>=4})

[1.65425211895987911, [c = HFloat(4.405119428776547)]]

> 

 


 

Download find_side_of_square_V2.mw

 

reply...

nm 9934
October 03 2024
0 0

Change your alpha to start from say 0.05 instead from 0 then it works.

You can see from your formula for M that you are dividing by alpha (inside the integral). Hence you can not use alpha=0 in the slider.
 

> 

params := {alpha = 1, gg = 0.1, k = 1, mu = 10, sigma = 5, w = 2};

{alpha = 1, gg = .1, k = 1, mu = 10, sigma = 5, w = 2}

> 

M := sqrt(-(gamma*(gamma*k*mu - 2*k*sigma)*k/(gamma*sigma - 1) + k^2)/alpha)*sinh(2*(gamma*k*(2*gamma*k*w + 2*k^2 + sigma^2)/(gamma*sigma - 1) - 2*k*sigma)*t/(gamma*sigma - 1))*exp(((2*gamma*k*w + 2*k^2 + sigma^2)*t/(gamma*sigma - 1) + sigma*x)*I)/(cosh(2*(gamma*k*(2*gamma*k*w + 2*k^2 + sigma^2)/(gamma*sigma - 1) - 2*k*sigma)*t/(gamma*sigma - 1)) - 1);

(-(gamma*(gamma*k*mu-2*k*sigma)*k/(gamma*sigma-1)+k^2)/alpha)^(1/2)*sinh(2*(gamma*k*(2*gamma*k*w+2*k^2+sigma^2)/(gamma*sigma-1)-2*k*sigma)*t/(gamma*sigma-1))*exp(((2*gamma*k*w+2*k^2+sigma^2)*t/(gamma*sigma-1)+sigma*x)*I)/(cosh(2*(gamma*k*(2*gamma*k*w+2*k^2+sigma^2)/(gamma*sigma-1)-2*k*sigma)*t/(gamma*sigma-1))-1)

> 

Explore(
     plot3d(abs(M), t = -5 .. 5, x = -10 .. 10,
        view = -10 .. 10, grid = [150, 150],
        color = [blue], style = surface, adaptmesh = false,
        size = [500, 500]),

     alpha = 0.05 .. 1.0,
     w = -5.0 .. 5.0,
     mu = -5.0 .. 5.0,
     sigma = -5.0 .. 5.0,
     k = -5.0 .. 5.0,

     placement = right);
> 

 


 

Download change_alpha.mw

 

 

reply...

nm 9934
September 25 2024
0 2

one way is to use allvalues

> 

eq1:=y=1/2*(x-5)^2-6;
eq2:=y=3*x-2;

y = (1/2)*(x-5)^2-6

y = 3*x-2

> 

sol:=solve([eq1,eq2],{x,y})

{x = RootOf(_Z^2-16*_Z+17), y = 3*RootOf(_Z^2-16*_Z+17)-2}

> 

allvalues(sol)

{x = 8-47^(1/2), y = 22-3*47^(1/2)}, {x = 8+47^(1/2), y = 22+3*47^(1/2)}

 

 

Download allvalues.mw

reply...

nm 9934
September 21 2024
4 0

First option

Use dsolve(ode_2,Lie) and now it can now solve the second ode

> 

restart;
> 

interface(version);

`Standard Worksheet Interface, Maple 2024.1, Windows 10, June 25 2024 Build ID 1835466`

> 

ode_2:=y(x) = (x - 1)*(x - 3)*diff(y(x), x) + exp((x*y(x)^2 - 4*x*y(x) + 3*x - y(x)^2 + 4*y(x) - 1)/(x - 3));

y(x) = (x-1)*(x-3)*(diff(y(x), x))+exp((x*y(x)^2-4*x*y(x)+3*x-y(x)^2+4*y(x)-1)/(x-3))

> 

sol_2:=dsolve(ode_2); #can't solve it

> 

restart; #make sure to do restart as Maple remembers the above result otherwise and will not try to solve it again!
> 

ode_2:=y(x) = (x - 1)*(x - 3)*diff(y(x), x) + exp((x*y(x)^2 - 4*x*y(x) + 3*x - y(x)^2 + 4*y(x) - 1)/(x - 3));
sol_2:=dsolve(ode_2,Lie); #Now it can

y(x) = (x-1)*(x-3)*(diff(y(x), x))+exp((x*y(x)^2-4*x*y(x)+3*x-y(x)^2+4*y(x)-1)/(x-3))

y(x) = (RootOf(-Intat(1/(exp(_a^2-1)-2), _a = _Z)*(x-3)^(1/2)+c__1*(x-3)^(1/2)+(x-1)^(1/2))*(x-3)^(1/2)+2*(x-1)^(1/2))/(x-1)^(1/2)

> 

 


 

Download option_1_symmetry.mw

Second option

My guess is that it did not solve it for the second ode due to complexity of the integral that results from solving the canonical ode in S(R) coordinates.

This below is a way to force it to  solve it. Using the infinitesimals generated from the first ode, which are verified on the second ode using symtest, generate the ODE in canonical coordinates. i,e, in S(R).  

Then use the transformation from S(R) to y(x) to obtain the solution in y(x). Maple provides function DEtools:-canoni for that.

Here is a small proc I use to do these sort of things. Now it solves the second ode.
 

> 

restart;
> 

interface(version);

`Standard Worksheet Interface, Maple 2024.1, Windows 10, June 25 2024 Build ID 1835466`

> 

ode_1:=y(x) = (x - 1)*(x - 3)*diff(y(x), x) + f((x*y(x)^2 - 4*x*y(x) + 3*x - y(x)^2 + 4*y(x) - 1)/(x - 3));
syms_1:=DEtools:-symgen(ode_1);

y(x) = (x-1)*(x-3)*(diff(y(x), x))+f((x*y(x)^2-4*x*y(x)+3*x-y(x)^2+4*y(x)-1)/(x-3))

[_xi = (x-1)^(1/2)*(x-3)^(3/2), _eta = (x-3)^(1/2)*(y-2)/(x-1)^(1/2)]

> 

ode_2:=y(x) = (x - 1)*(x - 3)*diff(y(x), x) + exp((x*y(x)^2 - 4*x*y(x) + 3*x - y(x)^2 + 4*y(x) - 1)/(x - 3));

y(x) = (x-1)*(x-3)*(diff(y(x), x))+exp((x*y(x)^2-4*x*y(x)+3*x-y(x)^2+4*y(x)-1)/(x-3))

> 

#First check that xi,eta satisfy second ode

> 

DEtools:-symtest(syms_1,ode_2, y(x))

0

> 

#solve it by hand now
get_ODE_in_canonical_v2:=proc(ode::`=`,y::symbol,x::symbol,S::symbol,R::symbol,xi,eta)
  local infinitesimals:=[rhs(xi),rhs(eta)];
  local tr,itr,ODE;
  tr  := DEtools:-canoni(infinitesimals,y(x),S(R));
  itr := op(1,[solve(tr,{x,y(x)})]);
  ODE:= PDEtools:-dchange(itr,ode,[R,S(R)],simplify);
  ODE:= op(solve(ODE,{diff(S(R),R)}));
  RETURN(simplify(ODE,symbolic),tr);
end proc:
> 

canonc_ode,tr:=get_ODE_in_canonical_v2(ode_2,y,x,S,R,syms_1[1],syms_1[2]);

diff(S(R), R) = -1/(-2+exp((R-1)*(R+1))), {R = (y(x)-2)*(x-1)^(1/2)/(x-3)^(1/2), S(R) = -(x-1)^(1/2)/(x-3)^(1/2)}

> 

sol_canonc:=dsolve(canonc_ode);

S(R) = Int(-1/(-2+exp((R-1)*(R+1))), R)+c__1

> 

eval(sol_canonc,tr); #convert back to y(x) coordinates

-(x-1)^(1/2)/(x-3)^(1/2) = Intat(-1/(-2+exp((_a-1)*(_a+1))), _a = (y(x)-2)*(x-1)^(1/2)/(x-3)^(1/2))+c__1

> 

sol:=PDEtools:-Solve(%,y(x)); #find explicit solution

y(x) = (RootOf(-Intat(1/(-2+exp((_a-1)*(_a+1))), _a = _Z)*(x-3)^(1/2)+c__1*(x-3)^(1/2)+(x-1)^(1/2))*(x-3)^(1/2)+2*(x-1)^(1/2))/(x-1)^(1/2)

> 

odetest(sol,ode_2); #verify

0

> 

 


 

Download force_solution_using_symmetry_sept_21_2024.mw

ps. you might want to submit bug report to Maplesoft on this. I think it should have solved it without the method Lie being given explicitly.

 

reply...

nm 9934
September 19 2024
3 2

 

 

> 

restart;
> 

mylist := [[sqrt(-5*x^2 - 5*x - 1) = -4*x - 1, {-1/3, -2/7}], [sqrt(-5*x^2 - 5*x - 1) = x + 1, {-1/2, -2/3}], [sqrt(-5*x^2 - 5*x - 1) = 4*x + 3, {-2/3, -5/7}], [sqrt(-5*x^2 - 5*x + 4) = -5*x + 3, {1/2, 1/3}], [sqrt(-5*x^2 - 4*x + 2) = -5*x + 2, {1/3, 1/5}], [sqrt(-5*x^2 - 4*x + 2) = -2*x + 1, {-1/3, 1/3}]]:
> 

toX:=s->latex(s,output=string):
s:="\\begin{enumerate}[label=\\arabic*)]\n":
for item in mylist do
    s:=cat(s,"\\item $",toX(item[1])," $\\hfill Answer: $",toX(item[2]),"$\n"):
od:
s:=cat(s,"\\end{enumerate}\n"):
    
> 

printf("%s",s)

\begin{enumerate}[label=\arabic*)]
\item $\sqrt{-5 x^{2}-5 x -1} = -4 x -1 $\hfill Answer: $\left\{-{\frac{2}{7}}, -{\frac{1}{3}}\right\}$
\item $\sqrt{-5 x^{2}-5 x -1} = x +1 $\hfill Answer: $\left\{-{\frac{2}{3}}, -{\frac{1}{2}}\right\}$
\item $\sqrt{-5 x^{2}-5 x -1} = 4 x +3 $\hfill Answer: $\left\{-{\frac{5}{7}}, -{\frac{2}{3}}\right\}$
\item $\sqrt{-5 x^{2}-5 x +4} = -5 x +3 $\hfill Answer: $\left\{{\frac{1}{2}}, {\frac{1}{3}}\right\}$
\item $\sqrt{-5 x^{2}-4 x +2} = -5 x +2 $\hfill Answer: $\left\{{\frac{1}{3}}, {\frac{1}{5}}\right\}$
\item $\sqrt{-5 x^{2}-4 x +2} = -2 x +1 $\hfill Answer: $\left\{-{\frac{1}{3}}, {\frac{1}{3}}\right\}$
\end{enumerate}

> 

 


 

Download to_latex_sept_19_2024.mw

reply...

nm 9934
September 17 2024
3 0

Basically, RootOf remembers something it calculated from first run which gets the error. In second run it uses that value it remmebered and this causes the error to go away.

You can see this by issuing trace(RootOf) and then in the second run you will see it says

       value remembered (at top level):

Worksheet below.

To make the error show each time, you can forget the remember table for RootOf. But I am not familar with remember tables so may be someone knows how to clear this specific table for RootOf. When I did  forget(RootOf) it had no effect and the error was still not showing up in the  second call.

May be I am not using the correct table in the forget command. But you get the idea. 

> 

restart;
> 

trace(RootOf);

RootOf

> 

RootOf(9*x1-5+RootOf(8*_Z^2+_Z-43, 41629632769253767815/18446744073709551616 .. 20814816384626883921/9223372036854775808), x1, 171590466306199/562949953421312 .. 343180932612401/1125899906842624);

{--> enter RootOf, args = 8*_Z^2+_Z-43, 41629632769253767815/18446744073709551616 .. 20814816384626883921/9223372036854775808

41629632769253767815/18446744073709551616, 20814816384626883921/9223372036854775808

false

8*_Z^2+_Z-43

{_Z}

{}

[]

8*_Z^2+_Z-43

8*z^2+z-43

FAIL

8*z^2+z-43

1

8*z^2+z-43

8*_Z^2+_Z-43

false

<-- exit RootOf (now at top level) = RootOf(8*_Z^2+_Z-43, 41629632769253767815/18446744073709551616 .. 20814816384626883921/9223372036854775808)}

{--> enter RootOf, args = 9*x1-5+RootOf(8*_Z^2+_Z-43, 41629632769253767815/18446744073709551616 .. 20814816384626883921/9223372036854775808), x1, 171590466306199/562949953421312 .. 343180932612401/1125899906842624

{x1}

{--> enter RootOf, args = 9*_Z-5+RootOf(8*_Z^2+_Z-43, 41629632769253767815/18446744073709551616 .. 20814816384626883921/9223372036854775808), 171590466306199/562949953421312 .. 343180932612401/1125899906842624

171590466306199/562949953421312, 343180932612401/1125899906842624

false

9*_Z-5+RootOf(8*_Z^2+_Z-43, 41629632769253767815/18446744073709551616 .. 20814816384626883921/9223372036854775808)

{_Z}

{RootOf(8*_Z^2+_Z-43, 41629632769253767815/18446744073709551616 .. 20814816384626883921/9223372036854775808)}

[RootOf(8*_Z^2+_Z-43, 41629632769253767815/18446744073709551616 .. 20814816384626883921/9223372036854775808) = rof[1]]

9*_Z-5+rof[1]

9*z-5+RootOf(8*_Z^2+_Z-43, 41629632769253767815/18446744073709551616 .. 20814816384626883921/9223372036854775808)

{--> enter RootOf, args = _Z^2-17

false

_Z^2-17

{_Z}

{}

[]

_Z^2-17

z^2-17

FAIL

z^2-17

1

z^2-17

_Z^2-17

false

<-- exit RootOf (now in \`evala/Indep/deg2\`) = RootOf(_Z^2-17)}

5/9-(1/9)*RootOf(8*_Z^2+_Z-43, 41629632769253767815/18446744073709551616 .. 20814816384626883921/9223372036854775808)

<-- ERROR in RootOf (now in RootOf) = cannot determine if this expression is true or false: %1, ln(0.1e11*abs(-.304805898398896+1.*RealRange(.304805898398895, .304805898398897)))/ln(10) < -6}

<-- ERROR in RootOf (now at top level) = cannot determine if this expression is true or false: %1, ln(0.1e11*abs(-.304805898398896+1.*RealRange(.304805898398895, .304805898398897)))/ln(10) < -6}

Error, (in property/ProbablyNonZero) cannot determine if this expression is true or false: ln(.1e11*abs(-.304805898398896+1.*RealRange(.304805898398895,.304805898398897)))/ln(10) < -6

> 

RootOf(9*x1-5+RootOf(8*_Z^2+_Z-43, 41629632769253767815/18446744073709551616 .. 20814816384626883921/9223372036854775808), x1, 171590466306199/562949953421312 .. 343180932612401/1125899906842624);

value remembered (at top level): RootOf(8*_Z^2+_Z-43, 41629632769253767815/18446744073709551616 .. 20814816384626883921/9223372036854775808) -> RootOf(8*_Z^2+_Z-43, 41629632769253767815/18446744073709551616 .. 20814816384626883921/9223372036854775808)

{--> enter RootOf, args = 9*x1-5+RootOf(8*_Z^2+_Z-43, 41629632769253767815/18446744073709551616 .. 20814816384626883921/9223372036854775808), x1, 171590466306199/562949953421312 .. 343180932612401/1125899906842624

{x1}

{--> enter RootOf, args = 9*_Z-5+RootOf(8*_Z^2+_Z-43, 41629632769253767815/18446744073709551616 .. 20814816384626883921/9223372036854775808), 171590466306199/562949953421312 .. 343180932612401/1125899906842624

171590466306199/562949953421312, 343180932612401/1125899906842624

false

9*_Z-5+RootOf(8*_Z^2+_Z-43, 41629632769253767815/18446744073709551616 .. 20814816384626883921/9223372036854775808)

{_Z}

{RootOf(8*_Z^2+_Z-43, 41629632769253767815/18446744073709551616 .. 20814816384626883921/9223372036854775808)}

[RootOf(8*_Z^2+_Z-43, 41629632769253767815/18446744073709551616 .. 20814816384626883921/9223372036854775808) = rof[1]]

9*_Z-5+rof[1]

9*z-5+RootOf(8*_Z^2+_Z-43, 41629632769253767815/18446744073709551616 .. 20814816384626883921/9223372036854775808)

5/9-(1/9)*RootOf(8*_Z^2+_Z-43, 41629632769253767815/18446744073709551616 .. 20814816384626883921/9223372036854775808)

<-- exit RootOf (now in RootOf) = 5/9-(1/9)*RootOf(8*_Z^2+_Z-43, 41629632769253767815/18446744073709551616 .. 20814816384626883921/9223372036854775808)}

<-- exit RootOf (now at top level) = 5/9-(1/9)*RootOf(8*_Z^2+_Z-43, 41629632769253767815/18446744073709551616 .. 20814816384626883921/9223372036854775808)}

5/9-(1/9)*RootOf(8*_Z^2+_Z-43, 41629632769253767815/18446744073709551616 .. 20814816384626883921/9223372036854775808)

 

 

Download twice_rootof.mw

reply...

nm 9934
September 16 2024
2 0

 

Maple is able to solve this analytically. Check if this solution is what you expected.
 

> 

restart;
> 

# the heat equation
pde := diff(v(x, t), t) = diff(v(x, t), x, x):

f := x -> exp(-x^2):

a := convert(-1.68858,rational):
b := convert(1.68858,rational):
IC:=v(a, t) = 0, v(b, t) = 0, v(x, 0) = 1 - int(f(f(s)), s = 0 .. x):
sol := pdsolve([pde,IC]);
sol2:=eval(rhs(sol),infinity=10):
F:=t0->evalf(eval(sol2,t=t0)):

v(x, t) = -(23669/39967)*(Sum(sin((1/79934)*n*Pi*(23669*x+39967))*exp(-(560221561/6389444356)*Pi^2*n^2*t)*(Int((-1+Int(exp(-exp(-2*s^2)), s = 0 .. x))*sin((1/79934)*n*Pi*(23669*x+39967)), x = -39967/23669 .. 39967/23669)), n = 1 .. infinity))

> 

plots[animate]( plot, [F(t), x=-5..5], t=0..2.5 );

> 

 


 

Download heat_pde.mw

reply...

nm 9934
September 16 2024
1 1

I do not think there is a way to tell dsolve to do this. But you can post-process the solution. Here is an example.

But pretty soon you will start having constants of integrations like c__100 and c__101 or may be c__5000 and so on depending on how many ODE's you plan to run in your loop. 
 

> 

interface(version);

`Standard Worksheet Interface, Maple 2024.1, Windows 10, June 25 2024 Build ID 1835466`

> 

restart;
> 

fix_my_C:=proc(sol::`=`)::`=`;
   local C_used::set,my_C::set;
   C_used:= indets(sol,And(symbol, suffixed(c__, nonnegint))):
   my_C := map(X->X=`tools/genglobal`(c__),C_used):
   eval(sol,my_C);
end proc:
 
> 

`tools/genglobal`[1](c__, 1, :-reset); #do once at start
> 

ode:=diff(y(x),x$3)+y(x)=sin(x);
sol:=dsolve(ode);
sol:=fix_my_C(sol);

diff(diff(diff(y(x), x), x), x)+y(x) = sin(x)

y(x) = -(1/2)*cos(x)/((3^(1/2)-2)*(2+3^(1/2)))-(1/2)*sin(x)/((3^(1/2)-2)*(2+3^(1/2)))+c__1*exp(-x)+c__2*exp((1/2)*x)*cos((1/2)*3^(1/2)*x)+c__3*exp((1/2)*x)*sin((1/2)*3^(1/2)*x)

y(x) = -(1/2)*cos(x)/((3^(1/2)-2)*(2+3^(1/2)))-(1/2)*sin(x)/((3^(1/2)-2)*(2+3^(1/2)))+c__1*exp(-x)+c__2*exp((1/2)*x)*cos((1/2)*3^(1/2)*x)+c__3*exp((1/2)*x)*sin((1/2)*3^(1/2)*x)

> 

ode:=diff(y(x),x$2)+y(x)=sin(x);
sol:=dsolve(ode);
sol:=fix_my_C(sol);

diff(diff(y(x), x), x)+y(x) = sin(x)

y(x) = sin(x)*c__2+cos(x)*c__1+(1/2)*sin(x)-(1/2)*cos(x)*x

y(x) = sin(x)*c__5+cos(x)*c__4+(1/2)*sin(x)-(1/2)*cos(x)*x

> 

ode:=diff(y(x),x)+y(x)=sin(x);
sol:=dsolve(ode);
sol:=fix_my_C(sol);

diff(y(x), x)+y(x) = sin(x)

y(x) = -(1/2)*cos(x)+(1/2)*sin(x)+c__1*exp(-x)

y(x) = -(1/2)*cos(x)+(1/2)*sin(x)+c__6*exp(-x)

> 

ode:=diff(y(x),x$5)+y(x)=0;
sol:=dsolve(ode);
sol:=fix_my_C(sol);

diff(diff(diff(diff(diff(y(x), x), x), x), x), x)+y(x) = 0

y(x) = c__1*exp(-x)-c__2*exp((-(1/4)*5^(1/2)+1/4)*x)*sin((1/4)*2^(1/2)*(5+5^(1/2))^(1/2)*x)-c__3*exp(((1/4)*5^(1/2)+1/4)*x)*sin((1/4)*2^(1/2)*(5-5^(1/2))^(1/2)*x)+c__4*exp((-(1/4)*5^(1/2)+1/4)*x)*cos((1/4)*2^(1/2)*(5+5^(1/2))^(1/2)*x)+c__5*exp(((1/4)*5^(1/2)+1/4)*x)*cos((1/4)*2^(1/2)*(5-5^(1/2))^(1/2)*x)

y(x) = c__7*exp(-x)-c__8*exp((-(1/4)*5^(1/2)+1/4)*x)*sin((1/4)*2^(1/2)*(5+5^(1/2))^(1/2)*x)-c__9*exp(((1/4)*5^(1/2)+1/4)*x)*sin((1/4)*2^(1/2)*(5-5^(1/2))^(1/2)*x)+c__10*exp((-(1/4)*5^(1/2)+1/4)*x)*cos((1/4)*2^(1/2)*(5+5^(1/2))^(1/2)*x)+c__11*exp(((1/4)*5^(1/2)+1/4)*x)*cos((1/4)*2^(1/2)*(5-5^(1/2))^(1/2)*x)

> 

ode:=diff(y(x),x$2)+y(x)=sin(x)+1;
sol:=dsolve(ode);
sol:=fix_my_C(sol);

diff(diff(y(x), x), x)+y(x) = sin(x)+1

y(x) = sin(x)*c__2+cos(x)*c__1+(1/2)*sin(x)+1-(1/2)*cos(x)*x

y(x) = sin(x)*c__13+cos(x)*c__12+(1/2)*sin(x)+1-(1/2)*cos(x)*x

> 

 


 

Download different_C_for_each_dsolve.mw

reply...

nm 9934
September 15 2024
3 1

If you convert it to poly then it works. The command is convert(p,polynom)

> 

p := asympt(x*1/(1 - a*x - b*x^2), x);
p:=convert(p,polynom);

-1/(b*x)+a/(b^2*x^2)-(1/b+a^2/b^2)/(b*x^3)-(-a/b^2-(a^2+b)*a/b^3)/(b*x^4)-((a^2+b)/b^3+a^2*(a^2+2*b)/b^4)/(b*x^5)+O(1/x^6)

-1/(b*x)+a/(b^2*x^2)-(1/b+a^2/b^2)/(b*x^3)-(-a/b^2-(a^2+b)*a/b^3)/(b*x^4)-((a^2+b)/b^3+a^2*(a^2+2*b)/b^4)/(b*x^5)

> 

coeff(p,1/x)

-1/b

 

 

Download coeffs.mw

reply...

nm 9934
September 14 2024
1 6

DESol in solution really means there is no closed form solution. At least Maple could not find one.

You did not give the input you used so can't verify.

But you can always use something like select or evalindets to pick these out. Below worksheet showing how. Maple 2024.1 on windows

> 

interface(version);

`Standard Worksheet Interface, Maple 2024.1, Windows 10, June 25 2024 Build ID 1835466`

> 

ode:=diff(y(x),x)-y(x)*a-b*y(x)^2=f(x);
sol:=dsolve(ode);

diff(y(x), x)-y(x)*a-b*y(x)^2 = f(x)

> 

my_sol := DESol( ode, y(x) );

DESol({diff(y(x), x)-y(x)*a-b*y(x)^2-f(x)}, {y(x)})

> 

#to extract just first argument of DESol
evalindets(my_sol,'specfunc'(anything,DESol),F->op(1,F))

{diff(y(x), x)-y(x)*a-b*y(x)^2-f(x)}

> 

#to extract all arguments
evalindets(my_sol,'specfunc'(anything,DESol),F->op(1..,F))

{diff(y(x), x)-y(x)*a-b*y(x)^2-f(x)}, {y(x)}

> 

 


 

Download pick_DESol_argument.mw

reply...

nm 9934
September 12 2024
3 3

There are two main issues. First from software design point of view, it is not good idea to couple two separate functionalities into one function. dsolve() should just solve an ode and nothing else. odetest() should just check that a solution verifies the ode.  

Even if it is optional for dsolve to test its solution before, I do not think it will be good to combine them. Too much complexity.  Also, dsolve in the process of solving an ode, could need to solve auxillary ode's.  Does it need to internally call odetest on each one of these?

But it is very easy to make your own my_dsolve() function which first calls dsolve, then calls odetest to verify maple solution.

This is what I do. Here is a basic algorithm how this can be done

my_dsolve:=proc(ode, IC, any other options.....)
  
    maple_solution:= dsolve( [ode,IC], any other options given ,..);
    if maple_solution not null then 
            #might need map here if more than one solution
            the_residual := map(X->odetest(X,[ode,IC]),[maple_solution]);
            if all the_residual entries are zero's then 
                 RETURN(maple_solution);
            else
                 RETURN(FAIL); #or return solution but with warning message. You decide.
            end if;
    else
           RETURN( maple_solution);  #nothing to verify. No solution
   end if;

end proc;

 

The second and more important issue is that odetest() is not always guaranteed to work.

It can hang on hard solutions, and so now dsolve have to deal with this new problem. How much timeoutput to use?  How does it know how much to wait? It also means dsolve will now take minutes to return a solution instead of seconds.

It also would mean it will not return solutions to ode's which it does now return a solution because it could not verify the solution (even though the solution could most likely be correct). I do not think people will be happy with this change.

odetest can also return false negative. i.e. it can return FAIL on solutions which is correct. see examples below. 

It is also possible that assumptions are need to verify the solution is correct or not. Sometimes I have to try 20 different assumptions to find one which makes odetest give zero.

But anything you do, do not use coulditbe() on result on odetest. This can easily give false positive. i.e. saying the residual can be zero, but the solution is wrong. But coulditbe() checks if the residual can be zero at just one point. 

Here are 10 examples showing the type of problems. I have hundreds more like this.

So making dsolve() do odetest() and return solution only if odetest verifies the solution is probably not something Maplesoft will want to do.

I think it is better that Maplesoft work more on improving odetest() itself and make it more robust but keep dsolve and odetest functions separate.

The problem of showing a solution verifies an ode is mathematically not an easy one. It is easy one for simple odes and simple solutions. But not so on much more complicated cases.

At school, the teacher says to just plug into the ode the solution, and then see if we get zero on both sides of the equation. But in real life this is not always easy to check if something is zero or not.

Worksheet attached
 

> 

 

Example 1

 
> 

ode:=(1+diff(y(x),x)^2)*(arctan(diff(y(x),x))+a*x)+diff(y(x),x) = 0;
sol:=dsolve(ode,y(x), singsol=all);

(1+(diff(y(x), x))^2)*(arctan(diff(y(x), x))+a*x)+diff(y(x), x) = 0

y(x) = Int(tan(RootOf(a*x*tan(_Z)^2+tan(_Z)^2*_Z+a*x+tan(_Z)+_Z)), x)+c__1

> 

odetest(sol,ode); #does not verify as is. Does this mean the solution is wrong or it just failed to verify?

(-arctan(tan(RootOf(2*a*x+sin(2*_Z)+2*_Z)))+RootOf(2*a*x+sin(2*_Z)+2*_Z))*tan(RootOf(2*a*x+sin(2*_Z)+2*_Z))/(a*x+RootOf(2*a*x+sin(2*_Z)+2*_Z))
> 

 

Example 2

 
> 

ode:=diff(y(x),x$2)^3=12*diff(y(x),x)*(x*diff(y(x),x$2)-2*diff(y(x),x));
sol:=dsolve(ode,y(x));

(diff(diff(y(x), x), x))^3 = 12*(diff(y(x), x))*(x*(diff(diff(y(x), x), x))-2*(diff(y(x), x)))

y(x) = (1/9)*x^4+c__1, y(x) = c__1, y(x) = Int(RootOf(-6*ln(x)-Intat((3*_f*(1/(_f*(9*_f-4)))^(1/2)*2^(1/3)*((3*(1/(_f*(9*_f-4)))^(1/2)*_f+1)^2*(9*_f-4)^4)^(1/3)-2*2^(2/3)*((3*(1/(_f*(9*_f-4)))^(1/2)*_f+1)*(9*_f-4)^2)^(1/3)-2^(1/3)*((3*(1/(_f*(9*_f-4)))^(1/2)*_f+1)^2*(9*_f-4)^4)^(1/3)+18*_f-8)/(_f*(9*_f-4)), _f = _Z)+6*c__1)*x^3, x)+c__2

> 

map(X->odetest(X,ode),[sol]): #does not verify 3rd solution and very slow to verify also. Too large to show residual

> 

 

Example 3

 
> 

ode:=(1+x^2)*diff(y(x),x$2)+1+diff(y(x),x)^2=x;
sol:=dsolve(ode,y(x)): #solution too large to show

(x^2+1)*(diff(diff(y(x), x), x))+1+(diff(y(x), x))^2 = x

> 

odetest(sol,ode);  #hangs or take too long to find out

> 

 

Example 4

 
> 

ode:=x^2*diff(y(x),x$2)+x^2*diff(y(x),x)+y(x)=0:
sol:=dsolve(ode,y(x),type='series',x=0):

> 

odetest(sol,ode,'series','point'=0);  #well known internal bug, been there for years. Fails to verify

Error, (in odetest/series) complex argument to max/min: 13/2-1/2*I*3^(1/2)

> 

 

> 

 

Example 5

 
> 

ode:=t^3*diff(y(t),t$2)+sin(t^3)*diff(y(t),t)+t*y(t)=0;
sol:=dsolve(ode,y(t),type='series',t=0):
odetest(sol,ode,'series','point'=0);  #well known internal bug, been there for years. Fails to verify

t^3*(diff(diff(y(t), t), t))+sin(t^3)*(diff(y(t), t))+t*y(t) = 0

Error, (in odetest/series) need to determine the sign of I*3^(1/2)

> 

 

> 

 

Example 6

 
> 

ode:=x^2*(1+x)*diff(y(x),x$2)-x*(1-6*x-x^2)*diff(y(x),x)+(1+6*x+x^2)*y(x)=0;
sol:=dsolve(ode,y(x),type='series',x=0);

x^2*(x+1)*(diff(diff(y(x), x), x))-x*(-x^2-6*x+1)*(diff(y(x), x))+(x^2+6*x+1)*y(x) = 0

y(x) = c__1*x*(series(1-12*x+(119/2)*x^2-(583/3)*x^3+(1981/4)*x^4-(80287/75)*x^5+O(x^6),x,6))+c__2*(x*ln(x)*(series(1-12*x+(119/2)*x^2-(583/3)*x^3+(1981/4)*x^4-(80287/75)*x^5+O(x^6),x,6))+x*(series(17*x-(471/4)*x^2+445*x^3-(118285/96)*x^4+(702451/250)*x^5+O(x^6),x,6)))

> 

odetest(sol,ode,'series','point'=0);  #fails to verify, but solution is correct

Warning, unable to compute series necessary to test the given solution

FAIL

> 

 

Example 7

 
> 

ode:=2*x^2*(2+x)*diff(y(x),x$2)+5*x^2*diff(y(x),x)+(1+x)*y(x)=0;
sol:=dsolve(ode,y(x),type='series',x=0);

2*x^2*(2+x)*(diff(diff(y(x), x), x))+5*x^2*(diff(y(x), x))+(x+1)*y(x) = 0

y(x) = c__1*x^(1/2)*(series(1-(3/4)*x+(15/32)*x^2-(35/128)*x^3+(315/2048)*x^4-(693/8192)*x^5+O(x^6),x,6))+c__2*(x^(1/2)*ln(x)*(series(1-(3/4)*x+(15/32)*x^2-(35/128)*x^3+(315/2048)*x^4-(693/8192)*x^5+O(x^6),x,6))+x^(1/2)*(series((1/4)*x-(13/64)*x^2+(101/768)*x^3-(641/8192)*x^4+(7303/163840)*x^5+O(x^6),x,6)))

> 

odetest(sol,ode,'series','point'=0);  #fails to verify, but solution is correct

Warning, unable to compute series necessary to test the given solution

FAIL

> 

 

> 

 

Example 8

 
> 

ode:=x^2*(2-x^2)*diff(y(x),x$2)-2*x*(1+2*x^2)*diff(y(x),x)+(2-2*x^2)*y(x)=0;
sol:=dsolve(ode,y(x),type='series',x=0);

x^2*(-x^2+2)*(diff(diff(y(x), x), x))-2*x*(2*x^2+1)*(diff(y(x), x))+(-2*x^2+2)*y(x) = 0

y(x) = c__1*x*(series(1+(3/4)*x^2+(15/32)*x^4+O(x^6),x,6))+c__2*(x*ln(x)*(series(1+(3/4)*x^2+(15/32)*x^4+O(x^6),x,6))+x*(series(-(1/8)*x^2-(13/128)*x^4+O(x^6),x,6)))

> 

odetest(sol,ode,'series','point'=0);  #fails to verify, but solution is correct

Warning, unable to compute series necessary to test the given solution

FAIL

> 

 

Example 9

 
> 

ode:=y(x) = arcsin(diff(y(x), x)) + ln(1 + diff(y(x), x)^2);
sol:=dsolve(ode,y(x));

y(x) = arcsin(diff(y(x), x))+ln(1+(diff(y(x), x))^2)

x-Intat(1/sin(RootOf(-_a+_Z+ln(sin(_Z)^2+1))), _a = y(x))-c__1 = 0

> 

odetest(sol,ode); #does not give zero, but solution is correct

-arcsin(sin(RootOf(-y(x)+_Z+ln(3/2-(1/2)*cos(2*_Z)))))+RootOf(-y(x)+_Z+ln(3/2-(1/2)*cos(2*_Z)))

> 

 

Example 10

 
> 

ode:=(1 + diff(y(x), x)^2)*(arctan(diff(y(x), x)) + a*x) + diff(y(x), x) = 0;
sol:=dsolve(ode,y(x));

(1+(diff(y(x), x))^2)*(arctan(diff(y(x), x))+a*x)+diff(y(x), x) = 0

y(x) = Int(tan(RootOf(a*x*tan(_Z)^2+tan(_Z)^2*_Z+a*x+tan(_Z)+_Z)), x)+c__1

> 

odetest(sol,ode); #does not give zero, but solution is correct

(-arctan(tan(RootOf(2*a*x+sin(2*_Z)+2*_Z)))+RootOf(2*a*x+sin(2*_Z)+2*_Z))*tan(RootOf(2*a*x+sin(2*_Z)+2*_Z))/(a*x+RootOf(2*a*x+sin(2*_Z)+2*_Z))

> 

 

> 

 

> 

 


 

Download on_odetest_sept_12_2024.mw

reply...

nm 9934
September 06 2024
4 4

Update

Here is an updated version wich now gets exact result in book. I used applyrule few places to force specific transformation.  Book uses C and D for constants. I used c__2 and c__3.

 

The idea is to replace G'/G^2  by u(t) and solve the resulting ode then replace u(t) back by G'/G^2

 

> 

restart;
> 

eq_47:= diff( diff(G(zeta),zeta)/G(zeta)^2,zeta)=mu+lambda*(diff(G(zeta),zeta)/G(zeta)^2)^2;
the_sub_1:=diff(G(zeta),zeta)/G(zeta)^2 = u(zeta);
the_sub_2:=dsolve(the_sub_1);
PDEtools:-dchange({the_sub_2},eq_47,{u},params={lambda,mu});
new_ode:=simplify(%);
u_sol:=dsolve(new_ode);
u_sol:=lhs(u_sol)=applyrule(tan(x::anything)=sin(x)/cos(x),rhs(u_sol));
 

(diff(diff(G(zeta), zeta), zeta))/G(zeta)^2-2*(diff(G(zeta), zeta))^2/G(zeta)^3 = mu+lambda*(diff(G(zeta), zeta))^2/G(zeta)^4

(diff(G(zeta), zeta))/G(zeta)^2 = u(zeta)

G(zeta) = 1/(Int(-u(zeta), zeta)+c__1)

(2*u(zeta)^2/(Int(-u(zeta), zeta)+c__1)^3+(diff(u(zeta), zeta))/(Int(-u(zeta), zeta)+c__1)^2)*(Int(-u(zeta), zeta)+c__1)^2-2*u(zeta)^2/(Int(-u(zeta), zeta)+c__1) = mu+lambda*u(zeta)^2

diff(u(zeta), zeta) = mu+lambda*u(zeta)^2

u(zeta) = tan((mu*lambda)^(1/2)*(c__1+zeta))*(mu*lambda)^(1/2)/lambda

u(zeta) = sin((mu*lambda)^(1/2)*(c__1+zeta))*(mu*lambda)^(1/2)/(cos((mu*lambda)^(1/2)*(c__1+zeta))*lambda)

> 

rules:=[ 'sin(x::anything)=sin('expand'(x)), cos(x::anything)=cos('expand'(x))'];
u_sol:=lhs(u_sol)=applyrule(rules,rhs(u_sol));

[sin(x::anything) = sin(('expand')(x)), cos(x::anything) = cos(('expand')(x))]

u(zeta) = sin(c__1*(mu*lambda)^(1/2)+zeta*(mu*lambda)^(1/2))*(mu*lambda)^(1/2)/(cos(c__1*(mu*lambda)^(1/2)+zeta*(mu*lambda)^(1/2))*lambda)

> 

rules:=[ sin(a::anything+b::anything)=sin(a)*cos(b)+sin(b)*cos(a), cos(a::anything+b::anything)=cos(a)*cos(b)-sin(a)*sin(b)]:
u_sol:=lhs(u_sol)=simplify(applyrule(rules,rhs(u_sol)));

u(zeta) = (sin(c__1*(mu*lambda)^(1/2))*cos(zeta*(mu*lambda)^(1/2))+sin(zeta*(mu*lambda)^(1/2))*cos(c__1*(mu*lambda)^(1/2)))*(mu*lambda)^(1/2)/((cos(c__1*(mu*lambda)^(1/2))*cos(zeta*(mu*lambda)^(1/2))-sin(c__1*(mu*lambda)^(1/2))*sin(zeta*(mu*lambda)^(1/2)))*lambda)

> 

u_sol:=eval(u_sol,[sin(c__1*sqrt(mu*lambda))=c__2,cos(c__1*sqrt(mu*lambda))=c__3]);
u_sol:=eval(u_sol,rhs(the_sub_1)=lhs(the_sub_1));
u_sol:=lhs(u_sol)=applyrule(sqrt(mu*lambda)/lambda=sqrt(mu/lambda),rhs(u_sol));

u(zeta) = (c__2*cos(zeta*(mu*lambda)^(1/2))+sin(zeta*(mu*lambda)^(1/2))*c__3)*(mu*lambda)^(1/2)/((c__3*cos(zeta*(mu*lambda)^(1/2))-c__2*sin(zeta*(mu*lambda)^(1/2)))*lambda)

(diff(G(zeta), zeta))/G(zeta)^2 = (c__2*cos(zeta*(mu*lambda)^(1/2))+sin(zeta*(mu*lambda)^(1/2))*c__3)*(mu*lambda)^(1/2)/((c__3*cos(zeta*(mu*lambda)^(1/2))-c__2*sin(zeta*(mu*lambda)^(1/2)))*lambda)

(diff(G(zeta), zeta))/G(zeta)^2 = (mu/lambda)^(1/2)*(c__2*cos(zeta*(mu*lambda)^(1/2))+sin(zeta*(mu*lambda)^(1/2))*c__3)/(c__3*cos(zeta*(mu*lambda)^(1/2))-c__2*sin(zeta*(mu*lambda)^(1/2)))

> 

 


 

Download simplification.mw

 

This is for Eq 49. Same thing but using hyperbolic trig relations

> 

restart;
> 

eq_47:= diff( diff(G(zeta),zeta)/G(zeta)^2,zeta)=mu+lambda*(diff(G(zeta),zeta)/G(zeta)^2)^2;
the_sub_1:=diff(G(zeta),zeta)/G(zeta)^2 = u(zeta);
the_sub_2:=dsolve(the_sub_1);
PDEtools:-dchange({the_sub_2},eq_47,{u},params={lambda,mu});
new_ode:=simplify(%);
u_sol:=dsolve(new_ode);
u_sol:=simplify(u_sol) assuming lambda*mu<0;

(diff(diff(G(zeta), zeta), zeta))/G(zeta)^2-2*(diff(G(zeta), zeta))^2/G(zeta)^3 = mu+lambda*(diff(G(zeta), zeta))^2/G(zeta)^4

(diff(G(zeta), zeta))/G(zeta)^2 = u(zeta)

G(zeta) = 1/(Int(-u(zeta), zeta)+c__1)

(2*u(zeta)^2/(Int(-u(zeta), zeta)+c__1)^3+(diff(u(zeta), zeta))/(Int(-u(zeta), zeta)+c__1)^2)*(Int(-u(zeta), zeta)+c__1)^2-2*u(zeta)^2/(Int(-u(zeta), zeta)+c__1) = mu+lambda*u(zeta)^2

diff(u(zeta), zeta) = mu+lambda*u(zeta)^2

u(zeta) = tan((mu*lambda)^(1/2)*(c__1+zeta))*(mu*lambda)^(1/2)/lambda

u(zeta) = -tanh((-mu*lambda)^(1/2)*(c__1+zeta))*(-mu*lambda)^(1/2)/lambda

> 

u_sol:=lhs(u_sol)=applyrule(tanh(x::anything)=sinh(x)/cosh(x),rhs(u_sol));

 

u(zeta) = -sinh((-mu*lambda)^(1/2)*(c__1+zeta))*(-mu*lambda)^(1/2)/(cosh((-mu*lambda)^(1/2)*(c__1+zeta))*lambda)

> 

rules:=[ 'sinh(x::anything)=sinh('expand'(x)), cosh(x::anything)=cosh('expand'(x))'];
u_sol:=lhs(u_sol)=applyrule(rules,rhs(u_sol));

[sinh(x::anything) = sinh(('expand')(x)), cosh(x::anything) = cosh(('expand')(x))]

u(zeta) = -sinh((-mu*lambda)^(1/2)*c__1+(-mu*lambda)^(1/2)*zeta)*(-mu*lambda)^(1/2)/(cosh((-mu*lambda)^(1/2)*c__1+(-mu*lambda)^(1/2)*zeta)*lambda)

> 

rules:=[ sinh(a::anything+b::anything)=sinh(a)*cosh(b)+cosh(a)*sinh(b), cosh(a::anything+b::anything)=cosh(a)*cosh(b)+sinh(a)*sinh(b)]:
u_sol:=lhs(u_sol)=simplify(applyrule(rules,rhs(u_sol)));

u(zeta) = -(sinh((-mu*lambda)^(1/2)*c__1)*cosh((-mu*lambda)^(1/2)*zeta)+cosh((-mu*lambda)^(1/2)*c__1)*sinh((-mu*lambda)^(1/2)*zeta))*(-mu*lambda)^(1/2)/((cosh((-mu*lambda)^(1/2)*c__1)*cosh((-mu*lambda)^(1/2)*zeta)+sinh((-mu*lambda)^(1/2)*c__1)*sinh((-mu*lambda)^(1/2)*zeta))*lambda)

> 

u_sol:=eval(u_sol,[sinh(sqrt(-mu*lambda)*c__1)=c__2,cosh(c__1*sqrt(-mu*lambda))=c__3]);
u_sol:=eval(u_sol,rhs(the_sub_1)=lhs(the_sub_1));
 

u(zeta) = -(c__2*cosh((-mu*lambda)^(1/2)*zeta)+c__3*sinh((-mu*lambda)^(1/2)*zeta))*(-mu*lambda)^(1/2)/((c__3*cosh((-mu*lambda)^(1/2)*zeta)+c__2*sinh((-mu*lambda)^(1/2)*zeta))*lambda)

(diff(G(zeta), zeta))/G(zeta)^2 = -(c__2*cosh((-mu*lambda)^(1/2)*zeta)+c__3*sinh((-mu*lambda)^(1/2)*zeta))*(-mu*lambda)^(1/2)/((c__3*cosh((-mu*lambda)^(1/2)*zeta)+c__2*sinh((-mu*lambda)^(1/2)*zeta))*lambda)

> 

 


 

Download simplification_rest.mw

 

For eq (50), I get zero in RHS. I do not know how book got what they show there.

 

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