Advanced Search

nm

11353 Reputation

20 Badges

13 years, 8 days
Contact nm

MaplePrimes Activity


These are answers submitted by nm

reply...

nm 11353
October 25 2024
4 13
> 

interface(version);

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

> 

data:=[[A,B],[0,B],[A,0],[0,0]]:
the_result :=Array(1..0):
ode:=diff(G(xi),xi)=G(xi)^2+A*G(xi)+B:
the_result ,= ["A","B","solution"]:

for item in data do
    eval(ode,[A=item[1],B=item[2]]):
    sol:=dsolve(%):
    the_result ,= [item[1],item[2],sol]:
od:
the_result:=convert(the_result,listlist):
DocumentTools:-Tabulate(the_result,weights=[10,10,100],width=30):
 

Download dsolve_oct_25_2024.mw

update

This version evaluates U=a0+a1*G(xi) after solving for G(xi) using current cases. And then verifies the U solution against the current updated F_ode.

May be this is what the question asks for?

> 

restart;
> 

interface(version);

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

> 

data:=[[A=A,B=B],[A=0,B=B],[A=A,B=0],[A=0,B=0]]:
the_result :=Array(1..0):
G_ode:=diff(G(xi),xi)=G(xi)^2+A*G(xi)+B:
F_ode := (-delta*eta^2 + alpha*eta)*diff(diff(U(xi), xi), xi) - U(xi)*(2*eta*gamma*theta*(delta*eta - alpha)*U(xi)^2 + eta^2*delta*k^2 + (-alpha*k^2 - 2*delta*k)*eta + 2*k*alpha + delta) = 0:
the_cases := [alpha = delta*(A^2*eta^2 - 2*eta^2*k^2 - 4*B*eta^2 + 4*eta*k - 2)/(A^2*eta - 2*eta*k^2 - 4*B*eta + 4*k),
             eta = eta,
             a[0] = -A/(2*gamma*theta*RootOf(_Z^2*gamma*theta + 1)),
             a[1] = RootOf(_Z^2*gamma*theta + 1)
             ]:
a0_original_value:=select(has,the_cases,a[0])[]:
a1_original_value:=select(has,the_cases,a[1])[]:
the_result ,= ["A","B",a0_original_value,a1_original_value,G(xi),U(xi)=a[0] + a[1]*G(xi),"odetest result"]:

for item in data do
    current_G_ode:=eval(G_ode,item):
    sol_G:=dsolve(current_G_ode):
    current_case := eval(the_cases,item):
    current_F_ode:=eval(F_ode,current_case):
    my_U_sol:= eval(a[0] + a[1]*G(xi),[op(current_case),sol_G]);
    was_verified:=odetest(U(xi)=my_U_sol,current_F_ode);
    a0_usesd :=rhs(select(has,current_case,a[0])[]);
    a1_usesd :=rhs(select(has,current_case,a[1])[]);
    the_result ,= [rhs(item[1]),rhs(item[2]),a0_usesd,a1_usesd,sol_G,my_U_sol,was_verified]:
od:
the_result:=convert(the_result,listlist):
DocumentTools:-Tabulate(the_result,weights=[10,10,50,50,100,100,20],width=80):
 


Download dsolve_oct_25_2024_V2.mw

reply...

nm 11353
October 20 2024
1 1

This is what I get in 2024.1.  Not sure if this is what you meant.

> 

M:=sin(x)/cos(x);
simplify(M);

sin(x)/cos(x)

tan(x)

> 

M:=1/sin(x);
simplify(M);

1/sin(x)

csc(x)

> 

M:=1/cos(x);
simplify(M);

1/cos(x)

sec(x)

> 

Q:=sqrt(beta/(B*cos(zeta*sqrt(-lambda))))

(beta/(B*cos(zeta*(-lambda)^(1/2))))^(1/2)

> 

simplify(Q) assuming lambda>0

(beta*sech(zeta*lambda^(1/2))/B)^(1/2)

 

 

Download trig_simplification.mw

reply...

nm 11353
October 15 2024
1 1

may be this is what you meant

> 

old_values := [alpha = 0.33101604, theta = -2.54098361, mu = 4.89071038, k = 5.0, A[1] = 2.70491803, a = 3.63387978];

[alpha = .33101604, theta = -2.54098361, mu = 4.89071038, k = 5.0, A[1] = 2.70491803, a = 3.63387978]

> 

new_values:=map(X->lhs(X)=MapleTA:-Builtin:-decimal(2,rhs(X)),old_values)

[alpha = .33, theta = -2.54, mu = 4.89, k = 5.00, A[1] = 2.70, a = 3.63]

> 

#or

> 

new_values:=map(X->lhs(X)=:-parse(sprintf("%.2f",rhs(X))),old_values)

[alpha = .33, theta = -2.54, mu = 4.89, k = 5.00, A[1] = 2.70, a = 3.63]

 

 

Download change_decimal.mw

Now use the new_values instead.

reply...

nm 11353
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 11353
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 11353
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 11353
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 11353
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 11353
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 11353
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 11353
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 11353
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 11353
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.

 

reply...

nm 11353
September 02 2024
1 5

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

 

reply...

nm 11353
August 28 2024
2 7

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 n in it to cancel. It only has sin(n*t) so you can't just cancel the n .

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

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