Preben Alsholm

13603 Reputation

22 Badges

19 years, 205 days

MaplePrimes Activity


These are answers submitted by Preben Alsholm

You have the code, so what's the problem:
 

restart;

ode:=diff(y(x),x)=1-2*y(x/2)^2;                                                                                                       
ics := y(0) = 2;

y3:= dsolve({ics, ode},numeric, delaymax = 1.7);

YY5:=plots:-odeplot(y3,x=0..10); 

This has been the case forever. Unassigned names are treated as nonzero by many Maple commands.
Take the simple examples:
 

restart;
solve(a*x=0,x); # 0
evalb(a=0); # false
is(a=0); # false
is(a=0) assuming a=0; # true, obviously
1/a;
1/0; # error

For solve there is these days the parametric option:
 

solve(a*x=0,x,parametric);

In your case solve/parametric could be used like this:

A:=Matrix([[1, 0, 0, 0, a__1], [0, 1, 0, 0, b__1], [-216, -18, 0, 0, c__1], [0, 0, 1, 0, d__1], [0, 0, 0, 1, e__1]]);
sys:=LinearAlgebra:-GenerateEquations(A,[x1,x2,x3,x4,x5],<0,0,0,0,0>);
solve(sys,[x1,x2,x3,x4,x5],parametric);

The answer shows correctly that the condition for the zero solution is:
1/216*c__1 + 1/12*b__1 + a__1 <> 0

I would begin by plotting your expression y which is clearly even, but imaginary for 1<x <2.
You don't need any readlib these days. It has been unnecessary for many years.

restart;

y:=sqrt((x^2*(x^2-4))/(x^2-1));
###
plot(y,x=0..10);
eval(y,x=3/2);
evalc(%);
limit(y,x=infinity);
limit(y,x=1,left);
ymax:=maximize(y,x=2..10);
ymin:=minimize(y,x=2..10);
ymin:=minimize(y,x=0..1-1e-10);

 

The expression  is singular at (x,y)=(0,0):
eval(R0,{x=0,y=0}); ## error
Try another point e.g. (x,y) = (1,0):
 

mtaylor(R0, [x=1, y=0], 5);

 

I think I managed to read your ode.
I would strongly advocate using odeplot instead of phaseportrait.
You cannot get a direction field (arrows) from odeplot, but your ode is non-autonomous, so a direction field has no meaning anyway.
 

restart;
ode:=diff(y(x),x,x)+8*diff(y(x),x)+9*y(x)+36*y(x)^2=2*sin(x);
## A plot of y(x) versus x:
DEtools:-phaseportrait(ode,y(x),x=-0.2..5,[[y(0)=1,D(y)(0)=0]],numpoints=1000);
# Rewriting ode as a first order system:
ode1:=diff(y(x),x)=v(x);
ode2:=subs(ode1,ode);
## A plot of diff(y(x),x) (i.e. v(x)) versus y(x):
DEtools:-phaseportrait([ode1,ode2],[y(x),v(x)],x=-.2..5,[[y(0)=1,v(0)=0]],numpoints=1000);
#### Much simpler to use dsolve/numeric directly on ode:
res:=dsolve({ode,y(0)=1,D(y)(0)=0},numeric);

## Now plot with odeplot:
plots:-odeplot(res,[y(x),diff(y(x),x)],-.2..5,numpoints=1000,thickness=2);

You may have a look at the help page ?evalf,details.
Also to get inspired you could try showstat(`evalf/constant/Catalan`).
Here is a rather silly example:

restart;
`evalf/constant/K`:=proc()  evalf(sqrt(2)*sin(sqrt(17))) end proc:
evalf(K);
evalf[15](K);
evalf[20](K);
evalf(exp(K)*Pi);
#################
showstat(`evalf/constant/Catalan`);

Since I believe that Maple's solve simply uses a formula for solving polynomials of degree less than 5 the answer that comes out from solve will just be the result of inputting the coefficients of the polynomial into the formula.
Using the very same polynomial as mmcdara, here is a graphical illustration, where the returned sequence from solve is kept in a list, not a set. Thereby we are not relying on any special knowledge of the ordering of sets in Maple.
In the grapical illustration you see that none of the 3 solutions remain real: The first in the list is red, the second blue, and the third is green.

restart;
pol:=x^3+a*x+1;
sol:=[solve](pol,x); # result: a list
with(plots):
complexplot(sol,a=-20..20,color=[red,blue,"DarkGreen"],thickness=3,style=point);
animate(complexplot,[sol,a=-20..A,color=[red,blue,"DarkGreen"],thickness=3,style=point], A=-20..20);

The animation follows   :

You seem to be assuming that (a^b)^c = a^(b*c).
But that isn't correct always: Take a = -1, b=2, and c=1/3:

((-1)^2)^(1/3); 
(-1)^(2/3); 
evalc(%);

a^b ( a <> 0) in complex analysis can be defined as exp(b*log(a)), where log is one of the branches of the logarithm. In Maple log is always chosen as the principal branch (and usually denoted by ln):

exp(2/3*ln(-1));

This isn't a procedure but it is simple:
 

d:=copy(c); # if you don't want to change c
d[..,1]:=floor~(d[..,1]);
d;
c;

If you want to change c:
 

c[..,1]:=floor~(c[..,1]);
c;

 

You get the animation using odeplot by adding frames=80:
 

p1:=odeplot(Solusix,[x(t),y(t)],0..20,numpoints=1000,color=green,frames=80):
p2:=odeplot(Solusi2,[x(t),y(t)],0..20,numpoints=1000,color=red,frames=80):
p3:=odeplot(Solusi3,[x(t),y(t)],0..20,numpoints=1000,color=blue,frames=80):
display(p1,p2,p3);

You have an error in your ode for n(t):
You have:
 

diff(n(t), t) = alpha_n*(1 - n(t)) - beta_n(t);

It should be:

diff(n(t), t) = alpha_n*(1 - n(t)) - beta_n;

After correcting for that just do:
 

dsolve(dsys1,numeric,method=classical[rk4]);

You can specify the stepsize as in:
 

res:=dsolve(dsys1,numeric,method=classical[rk4],stepsize=0.01);
res(2); #The result at t = 2

 

Correcting some syntactical problems I got this ODE:
 

ODE:= (2/3)*int(diff(y(x),x)*x^2/(x^2 -1),x) =int(-sqrt(2*y(x)),x);

Is this the one you intended?
If so then the answer in Maple 2023.1 (and also in Maple 2021.2 and Maple 12)  using dsolve is
 

sqrt(y(x)) + 3/4*sqrt(2)/x + 3/4*sqrt(2)*x - c__1 = 0

Differentiating ODE you get a separable ode:
 

ode2:=diff(ODE,x);

From which the resulting implicit solution follows by separation of variables and then integrating both sides wrt. x.

But my interpretation of your original ODE may be wrong?
 

To get the floats do:
 

indets(X1,float);

 

Save the plot itself in some variable. Below I use p.

p:=%; # saving the plot right after executing the plot3d command.
data:=plottools:-getdata(p);
data[1];
data[2];
data[3];
dx:=5./48; #spacing between x-values
dy:=evalf(b)/48; ##spacing between y-values
X:=[seq(dx*i,i=0..48)]; # List of x-values
Y:=[seq(dy*i,i=0..48)]; # List of y-values
[X[7],Y[9]]; # An example xy point:
data[3][7,9]; # The corresponding value of phi[2]
eval(phi[2],{x=X[7],y=Y[9]}); # Check

As you notice this series is telescoping. So this works:
 

restart;
A:=sum(1/f(n)-1/f(n+1),n=1..infinity); # Using the unassigned f instead of sqrt
B:=subs(infinity=N,A);
expand(B);
eval(%,f=sqrt);
limit(%,N=infinity);

 

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