nm

11483 Reputation

20 Badges

13 years, 82 days

MaplePrimes Activity


These are answers submitted by nm

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

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.

 

I do not wish to find \n in the answer

    \n is special. It does not require escaping. So just change all your \\n  in your input to \n and you will not see \n anymore but will get new line instead, which is I assume what you wanted.

by adding extra \  you basically have scaped the next \. You do not want to do this with \n

update

Looks like OP just want to see "X","Y","buff"  each on its own row.

So why not just a Matrix? Forget about using printf for making display.


 

restart;

intersections := proc(P, Q, T)
local R, W, w, t, a, b, sol, buff, v;
local my_collection:=Array(1..0):
my_collection ,= "X","Y","buff";

sol := NULL;

if T = Y then
   W := X;
else
   W := Y;
fi;

R := resultant(P, Q, T);

#print(`Résultant :`); print(R);

w := fsolve(R, W);
t := NULL;

for v in [w] do
    t := t, fsolve(subs(W = v, P), T);
od;

for a in {w} do
    for b in {t} do
        if T = Y then
           buff := abs(subs(X = a, Y = b, P)) + abs(subs(X = a, Y = b, Q));
           my_collection ,= a,b,buff;

           #printf(`X=%a,   Y=%a   --->  %a\\n`, a, b, buff);
           if buff < 1/100000000 then
              sol := sol, [a, b];
           fi;
        else
            buff := abs(subs(X = b, Y = a, P)) + abs(subs(X = b, Y = a, Q));
            my_collection ,= a,b,buff;
            #printf(`X=%a,   Y=%a   --->  %a\\n `, a, b, buff);
            if buff < 1/100000000 then
               sol := sol, [b, a];
            fi;
        fi;
    od;
od;

#printf(`Nombre de solutions :  %a\\n`, nops({sol}));
#print({sol});

my_collection:=convert(my_collection,list);
my_collection:=Matrix( (nops(my_collection)/3,3),my_collection);
RETURN(my_collection);

end proc:

intersections(X^2 + Y^2 - 1, X - Y, X);

Matrix(5, 3, {(1, 1) = "X", (1, 2) = "Y", (1, 3) = "buff", (2, 1) = -.7071067812, (2, 2) = -.7071067812, (2, 3) = 0., (3, 1) = -.7071067812, (3, 2) = .7071067812, (3, 3) = 1.414213562, (4, 1) = .7071067812, (4, 2) = -.7071067812, (4, 3) = 1.414213562, (5, 1) = .7071067812, (5, 2) = .7071067812, (5, 3) = 0.})

 


 

Download display_stuff.mw

 

I do not see how the book solution your wrote can be the same as Maple's solution. 

In Maple solution, there is "n" in the denominator inside the sum. In your book solution you have removed this "n".

Now the "n" in Maple solution does cancel with one of the terms in the sum in the numerator, the one with cos(nt)*n but the second term in the numerator has no in it to cancel. It only has sin(n*t) so you can't just cancel the .

So your book solution is either wrong or you copied it wrong.

You can see below Maple solution and the book solution give different answer.  You can change infinity to as many terms your want and pick arbitrary x and t values and arbitray functions f(x) and g(x).

Btw, why do you not also insert content when you post your worksheet, instead of just adding link? It is small. 
 

 

 

restart;

 

# Definieer de parameters
L := Pi:
c := 1: # Je kunt c aanpassen aan je specifieke probleem

# Stel de golfvergelijking op
pde := diff(u(x,t),t,t) = c^2 * diff(u(x,t),x,x):

# Definieer de randvoorwaarden
bc := u(0,t) = 0, u(L,t) = 0:

# Definieer de beginvoorwaarden
ic := u(x,0) = f(x), D[2](u)(x,0) = g(x):

# Los de PDE op met begin- en randvoorwaarden
sol := pdsolve([pde, bc, ic], u(x,t)):

# Toon de oplossing
simplify(sol);

u(x, t) = 2*(Sum(sin(n*x)*((Int(sin(n*x)*f(x), x = 0 .. Pi))*cos(n*t)*n+(Int(sin(n*x)*g(x), x = 0 .. Pi))*sin(n*t))/n, n = 1 .. infinity))/Pi

eval(rhs(sol),infinity=5);
simplify(value(%)):
value(eval(%,[t=2,f(x)=x^3,g(x)=x^2])):
evalf(eval(%,x=1));

Sum(2*sin(n*x)*((Int(sin(n*x)*f(x), x = 0 .. Pi))*cos(n*t)*n+(Int(sin(n*x)*g(x), x = 0 .. Pi))*sin(n*t))/(Pi*n), n = 1 .. 5)

10.00905987

book_sol:=2*Sum(sin(n*x)*(Int(sin(n*x)*f(x), x = 0 .. Pi)*cos(n*t) + Int(sin(n*x)*g(x), x = 0 .. Pi)*sin(n*t)), n = 1 .. infinity)/Pi:
eval(book_sol,infinity=5);
simplify(value(%)):
value(eval(%,[t=2,f(x)=x^3,g(x)=x^2])):
evalf(eval(%,x=1));

2*(Sum(sin(n*x)*((Int(sin(n*x)*f(x), x = 0 .. Pi))*cos(n*t)+(Int(sin(n*x)*g(x), x = 0 .. Pi))*sin(n*t)), n = 1 .. 5))/Pi

12.43548172

 


 

Download not_same.mw

Use StringTools and make sure y is "y"

fx:=3*x^2-9*y;
str:=String(fx);
hasy:=StringTools:-Has(str,"y")

   true

The `->` confuses the translator. It is a `Rule` in Mathematica.  For correct translation use what Mathematica uses in its internal FullForm of the expression instead of the displayed Form. i.e replace `->` by `Rule` which is what it is.

So your code now works and becomes

restart

MmaTranslator:-FromMma("Subscript[x, 0] -> 2 (y - z)^2"); # should return: x[0] = 2*(y - z)^2

x[0] = (2*y-2*z)^2

MmaTranslator:-FromMma("Rule[ Subscript[x, 0],2 (y - z)^2]"); # should return: x[0] = 2*(y - z)^2

x[0] = 2*(y-z)^2

MmaTranslator:-FromMma("a -> 2*(y - z)^2");

a = (2*y-2*z)^2

MmaTranslator:-FromMma("Rule[a,2*(y - z)^2]");

a = 2*(y-z)^2

 


 

Download use_rule.mw

On a side note, I would avoid translating code with Subscript to Maple. Maple translates Subscript[x, 0] to indexed variable. Which on the screen looks like x__0 in Maple. But it is not the same. This can cause errors in your translated code if you mix it with your own Maple code which uses x__0 and not x[0].

i.e. Subscript[x, 0]  is translated to x[0] in Maple and not to x__0, and these are different types in Maple.

Here is example

restart;

e:=MmaTranslator:-FromMma(`Subscript[x, 0]`);

x[0]

whattype(e);

indexed

lprint(e);

x[0]

e:=x__0;

x__0

whattype(e);

symbol

lprint(e);

x__0


btw, I never ever use Subscripted variables in Mathematica. This has nothing to do with translation to Maple, but these cause their own share of problems in Mathematica itself. stackexchange Mathematica forum has dozens of questions related to using Subscripted variables in Mathematica and almost eveyone says to avoid them in Mathematica.

 

Download use_rule_V2.mw

 

plots:-polarplot(4, coordinateview = [0 .. 6,0 .. 2*Pi], thickness = 3,tickmarks=[default, 2*Pi/(Pi/12)]);

 

 

pde:=diff(u(x,t),x$2)-(v-2)*u(x,t)-(v-v1/mu)*(-u(x,t)^3+u(x,t)^5)=0

diff(diff(u(x, t), x), x)-(v-2)*u(x, t)-(v-v1/mu)*(-u(x, t)^3+u(x, t)^5) = 0

xx:=solve(zeta=sqrt(k)*(x-c*t),x);

(k^(1/2)*c*t+zeta)/k^(1/2)

tr:={x=xx,u(x,t)=U(zeta)};

 

{x = (k^(1/2)*c*t+zeta)/k^(1/2), u(x, t) = U(zeta)}

pde2:=PDEtools:-dchange(tr,pde,{U,zeta},'known'={u(x,t)},'unknown'={U(zeta)})

(diff(diff(U(zeta), zeta), zeta))*k-(v-2)*U(zeta)-(v-v1/mu)*(-U(zeta)^3+U(zeta)^5) = 0

pde3:=PDEtools:-dchange({U(zeta)=sqrt(W(zeta))},pde2,{W},'known'={U(zeta)},'unknown'={W(zeta)})

(-(1/4)*(diff(W(zeta), zeta))^2/W(zeta)^(3/2)+(1/2)*(diff(diff(W(zeta), zeta), zeta))/W(zeta)^(1/2))*k-(v-2)*W(zeta)^(1/2)-(v-v1/mu)*(-W(zeta)^(3/2)+W(zeta)^(5/2)) = 0

expand(pde3*sqrt(W(zeta)))

-(1/4)*k*(diff(W(zeta), zeta))^2/W(zeta)+(1/2)*(diff(diff(W(zeta), zeta), zeta))*k-W(zeta)*v+2*W(zeta)+W(zeta)^2*v-W(zeta)^3*v-W(zeta)^2*v1/mu+W(zeta)^3*v1/mu = 0

collect(%,W(zeta));

(-v+v1/mu)*W(zeta)^3+(v-v1/mu)*W(zeta)^2+(-v+2)*W(zeta)+(1/2)*(diff(diff(W(zeta), zeta), zeta))*k-(1/4)*k*(diff(W(zeta), zeta))^2/W(zeta) = 0

 


 

Download transformation.mw

convert the result to set, this will automatically remove duplicates.


 

restart;

interface(rtablesize=70);


result:=Array(1..0):
phi := (p__1, p__2, p__3, p__4, q__1, q__2, q__3, q__4, xi) -> (p__1*cosh(xi) - p__2*sinh(xi))/(p__3*cosh(xi) + p__4*sinh(xi));
xi := x:
for p__1 in [1, -1, I, -I] do
    for p__2 in [1, -1, I, -I] do
        for p__3 in [1, -1, -I, I] do
            for p__4 in [0] do
                for q__1 in [1, -1, I, -I] do
                    for q__2 in [1, -1, I, -I] do
                        for q__3 in [1, -1, I, -I] do
                            for q__4 in [1, -1, I, -I] do
                                result ,=evalf(phi(p__1, p__2, p__3, p__4, q__1, q__2, q__3, q__4, xi)):
                                #print('phi'(p__1, p__2, p__3, p__4, q__1, q__2, q__3, q__4, xi) = result1);
                            od:
                       od:
                    od:
                od:
            od:
         od:
     od:
od:

[10, 10]

proc (p__1, p__2, p__3, p__4, q__1, q__2, q__3, q__4, xi) options operator, arrow; (p__1*cosh(xi)-p__2*sinh(xi))/(p__3*cosh(xi)+p__4*sinh(xi)) end proc

print("Before, the size was ",numelems(result));
result_no_duplidates:=convert(convert(result,set),list):
print("After remvoing duplicates, the size is ",nops(result_no_duplidates));
map(X->[''phi(p__1, p__2, p__3, p__4, q__1, q__2, q__3, q__4, xi)'' = X],result_no_duplidates):
Matrix(%);

"Before, the size was ", 16384

"After remvoing duplicates, the size is ", 64

Matrix(%id = 36893489808398873292)

 


 

Download A.mw

I do not think there is conversion builtin in Maple for this, at least I could not find one.  But see if this does what you want for now. I only tested it on your input !

e0:=int(x^(1/2),x);
evalindets(e0,sqrt,F->sqrt(op(1,F)^numer(op(2,F))))

Update

For more general cases, (for your other comment below also) you can try this function called mySimplify


 

restart;

mySimplify :=e->evalindets(e,`^`,F->`if`( op(2,F)::fraction, surd( op(1,F)^numer(op(2,F)), denom(op(2,F))),F)):
test_cases:=[  int( (x^2)^(1/3),x), int(x^(1/2),x), x^(2/3) ]:
Matrix(map(X->[X,mySimplify(X)], test_cases))

Matrix(3, 2, {(1, 1) = (3/5)*x*(x^2)^(1/3), (1, 2) = (3/5)*x*(x^2)^(1/3), (2, 1) = (2/3)*x^(3/2), (2, 2) = (2/3)*sqrt(x^3), (3, 1) = x^(2/3), (3, 2) = (x^2)^(1/3)})

 


 

Download mysimplify.mw

while running my code i am not able to edit it

If you mean you are running something in the worksheet, and while it is sitll running you can't edit or do anything? If so, then this is by design. Maple interface can't handle making any changes to worksheet while it is running something. You have to wait until computation is done.

Make sure you set each worksheet to separate engine. This is a nice feature in Maple.

This way while one sheet is busy, you can always open new worksheet and use that because each will be using new engine. See options on how to change this setting.

 

although my code is working for one pde but while changing pde it is giving me error why?

Good chance you have made some error typing the new pde or the new command. The error message you get should make it clear why.

 some times it takes alot of time to evaluate is there any error in the code or should i change my laptop?

You should use the command timelimit look it up in help. You set a timeout how long some computation should take. The general format for doing this can be something like

try
   result := timelimit(30, some_maple_command); #30 seconds timeout
catch:
  print("Warning, The command has timed out");
  result:= FAIL;
end try;

You do not need to get new laptop. Many CAS computations take long time or hang. So you should always use timelimit to guard against this. Ofcourse having faster PC with more memory is not going to hurt if you can afford it. 

 

PDEtools:-Infinitesimals

compute the infinitesimals of symmetry groups admitted by a given PDE system

25944

pde:= diff(u(t,x,y),t$2)+2*a(t)*diff(u(t,x,y),x)^2+2*b(t)*u(t,x,y)*diff(u(t,x,y),x$2)+diff(u(t,x,y),y$2)=0;
PDEtools:-Infinitesimals(pde,u(t,x,y))

diff(diff(u(t, x, y), t), t)+2*a(t)*(diff(u(t, x, y), x))^2+2*b(t)*u(t, x, y)*(diff(diff(u(t, x, y), x), x))+diff(diff(u(t, x, y), y), y) = 0

[_xi[t](t, x, y, u) = 0, _xi[x](t, x, y, u) = 1, _xi[y](t, x, y, u) = 0, _eta[u](t, x, y, u) = 0], [_xi[t](t, x, y, u) = 0, _xi[x](t, x, y, u) = 0, _xi[y](t, x, y, u) = 1, _eta[u](t, x, y, u) = 0], [_xi[t](t, x, y, u) = 0, _xi[x](t, x, y, u) = x, _xi[y](t, x, y, u) = 0, _eta[u](t, x, y, u) = 2*u]

 

 

Download lie_pde.mw

Do not know if this will work for what you are doing or not.
 

26000

restart;

26000

interface(typesetting=extended);

extended

make_nice:=proc(k::anything)
     evalindets(k,`&*`(anything,'specfunc(ln)'),F->InertForm:-MakeInert(F));
end proc;

proc (k::anything) evalindets(k, `&*`(anything, 'specfunc(ln)'), proc (F) options operator, arrow; InertForm:-MakeInert(F) end proc) end proc

k:=3/8*ln(55/52);

(3/8)*ln(55/52)

make_nice(k)

`%*`(`%/`(3, 8), %ln(`%/`(55, 52)))

k:=3/8*ln(55/52)+sin(x)+3/4*exp(x);

(3/8)*ln(55/52)+sin(x)+(3/4)*exp(x)

make_nice(k)

`%*`(`%/`(3, 8), %ln(`%/`(55, 52)))+sin(x)+(3/4)*exp(x)

 


 

Download make_nice.mw

From help it says subs does not work on contravariant indecies. Only on covariant. So need to use SubstituteTensorIndices from Physics package, See help page on SubstituteTensorIndices for more info.  It says to use subs on contra need to use ~=~ format then it should work.

I am using V 2024. If this does not work for you, may be it is version issue then,

20936

restart;

20936

with(Physics):
Setup(mathematicalnotation = true);
Define(A, B):
e:=g_[alpha, mu] * A[~mu] * g_[~alpha, ~nu] * B[nu, sigma, ~rho];

[mathematicalnotation = true]

`Defined objects with tensor properties`

Physics:-g_[alpha, mu]*A[`~mu`]*Physics:-g_[`~alpha`, `~nu`]*B[nu, sigma, `~rho`]

subs(alpha = beta,e); #does not change contra

Physics:-g_[beta, mu]*A[`~mu`]*Physics:-g_[`~alpha`, `~nu`]*B[nu, sigma, `~rho`]

SubstituteTensorIndices(alpha = beta,e); #changes contra and covariant

Physics:-g_[beta, mu]*A[`~mu`]*Physics:-g_[`~beta`, `~nu`]*B[nu, sigma, `~rho`]

#help says this should work.
subs(~alpha = ~beta,e); #changes contra only OK

Physics:-g_[alpha, mu]*A[`~mu`]*Physics:-g_[`~beta`, `~nu`]*B[nu, sigma, `~rho`]

 

 

Download tensor_june_14_2024.mw

I never heard of the method of isodine. May be you meant isoclines  ?

To make slope plot do

ode:= x+y(x)*diff(y(x),x)=0;
DEtools:-DEplot(ode,y(x),x=-2..2,y=-3..3)

To get the solution, the command is 

dsolve(ode)

To get step by step solution the command is

Student:-ODEs:-ODESteps(ode);

3 4 5 6 7 8 9 Last Page 5 of 20