Preben Alsholm

13663 Reputation

22 Badges

19 years, 344 days

MaplePrimes Activity


These are answers submitted by Preben Alsholm

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);

 

I shall use Maple 12 (since I don't have your Maple 13; later versions have the almost same example in the help for Statistics:-Fit).
B is the independent variable, Bdata is a Vector of values for B, Hdata is a Vector  of values for the dependent variable.

restart;
with(Statistics):
Bdata := Vector([1, 2, 3, 4, 5, 6], datatype = float);
Hdata := Vector([2, 3, 4.8, 10.2, 15.6, 30.9], datatype = float);
H:=Fit(a*B+b*sinh(c*B), Bdata, Hdata, B);

Result: H := 1.28528203195904544*B+.244370171891131000*sinh(.873796260764654664*B)

Plotting H versus B and plotting the data points:
 

plot(H,B=1..6,labels=[B,'H']); p1:=%:
plot(Bdata,Hdata,style=point,symbol=solidcircle,symbolsize=20,color=blue); p2:=%:
plots:-display(p1,p2);

Choosing as initial R, not 0, but here 0.1:
 

restart;
#with(plots):with(plottools):with(DETools):

Sys:=diff(T(R),R)=((1-1/R)*(sqrt(1-(alpha/R)^2*(1-1/R))))^(-1),diff(Phi(R),R)=(alpha/R)^2*(sqrt(1-(alpha/R)^2*(1-1/R)))^(-1);

inits:=[[T(0)=0.5,Phi(0)=0],[T(0)=0.5,Phi(0)=Pi/4]];
K:=dsolve([Sys,op(op(1,inits))],[Phi(R),T(R)],numeric,parameters=[alpha],output=listprocedure);
eval(Sys,R=0);
## Now changing from 0 to 0.1 in inits:
inits2:=[[T(0.1)=0.5,Phi(0.1)=0],[T(0.1)=0.5,Phi(.1)=Pi/4]];
K:=dsolve([Sys,op(op(1,inits2))],[Phi(R),T(R)],numeric,parameters=[alpha],output=listprocedure);
K(parameters=[1]); # Picking alpha = 1 (arbitrarily)
plots:-odeplot(K,[[R,T(R)],[R,Phi(R)]],0..1);

Notice that simplification helps some. It shows that only Phi has a problem at 0, but both have a singularity at 1:
 

simplify([Sys]) assuming R>0, alpha::real;

Output:

Proceeding with Sys2 and changing 0 to 1e-8 in the initialconditions:
 

Digits:=15:
inits3:=[[T(1e-8)=0.5,Phi(1e-8)=0],[T(1e-8)=0.5,Phi(1e-8)=Pi/4]];
K3:=dsolve([op(Sys2),op(op(1,inits3))],[Phi(R),T(R)],numeric,parameters=[alpha],output=listprocedure,abserr=1e-15,relerr=1e-10)
K3(parameters=[1]);
plots:-odeplot(K3,[[R,T(R)],[R,Phi(R)]],0..1);

 

Here is a version, where A is created as a table:
 

restart;
L:=[seq(2*i - 1, i = 1 .. 10)];
for i in L do
  A[i]:=(x/a)^(i + 1)*(1 - x/a)^2
end do;
eval(A);

 

MultiSeries:-limit(y(0),y(0)=1);
MultiSeries:-limit(y(x),y(x)=1);

These give a more reasonable error:

Error, invalid input: MultiSeries:-limit expects its 2nd argument, eq, to be of type name = algebraic, but received y(0) = 1

type(y(0),name); #  false

You could use frontend:
 

frontend(limit,[y(0),y(0)=1],[{`+`,`*`,`=`}]);

 

Use subsindets (not evalindets):
 

restart;
p:=piecewise(i <= 0., [], i < 0.7900000000, [{t2 = 2.5*RootOf((270343941*exp((25*_Z)/4)*i)/400 + (270343941*exp((75*_Z)/4)*i)/400 - (270343941*exp((25*_Z)/2)*i)/200 - (531032741*i*exp((29*_Z)/4))/500 + (347585067*i*exp((27*_Z)/2))/200 - (12873521*i*exp((79*_Z)/4))/16 + (16*exp((29*_Z)/4))/25 + (12873521*exp(_Z)*i)/100)}], 0.7900000000 <= i, []);
p1:=subsindets(p,list,s->op(s));

My p is your s2.

This following observation may not explain what you are seeing, but here it is:
 

restart;
whattype(ECS_8.3); #`*` (product)
lprint(ECS_8.3);
T[ECS_8.3]:=dot3;
eval(T);
indices(T);
lprint(%);
T[`ECS_8.3`]:=Name;
indices(T);
T["ECS_8.3"]:=maplestring;
eval(T);
indices(T,nolist,pairs);

The crucial point here is that   a.3 is taken to mean a*3 when a is a name, which ECS_8 is.

I think a better way is this:
 

restart;
int(1/(sqrt(x__0 - x)*sqrt(-x^2 + 1)), x = 0 .. x__0);
# So we make assumptions:
A:=int(1/(sqrt(x__0 - x)*sqrt(-x^2 + 1)), x = 0 .. x__0) assuming x__0>=0,x__0<=1;
f:=unapply(A,x__0);
plot(f(x__0), x__0 = 0..1 ,labels=[x__0,'f(x__0)']);
f(0.5);
h := y -> fsolve(f(x__0) = y, x__0 =0..1); 
h(f(0.5));

This was done in Maple 2023.0, but works in Maple 2021.2 and Maple 2022.2 too.

Since you can look up all the definitions required, it would be helpful if you could point out which result is wrong.

restart;
L:=[sqrt(-3)^1, sqrt(-3)^7]; # Output L := [sqrt(3)*I, -27*I*sqrt(3)]
?type/radnumext
eval(`type/radnumext`);
?type/radext
interface(verboseproc=3);
eval(`type/radext`);
?type/radical
?radnum
eval(`type/radnum`);

Continuing with L:

L := [sqrt(3)*I, -27*I*sqrt(3)]
LL:=convert~(L,list);
###
seq(type~(LL[i],radnumext),i=1..2);
map(indets, L, radnumext);
###
seq(type~(LL[i],radext),i=1..2);
map(indets, L, radext);
###
seq(type~(LL[i],radical),i=1..2);
map(indets, L, radical);

One might be surprised that converting -27*I*sqrt(3) to a list gives only 2 elements.
The reason is probably the definition of I as Complex(1). See ? I.
 

L[2];
nops(%);
op(1,L[2]);
op(2,L[2]);
Complex(-27);
addressof(op(1,L[2]));
addressof(Complex(-27));## The same address
 

Continuing a little further:
 

addressof(L[2]);
addressof(Complex(-27*sqrt(3))); # Not the same as L[2]
simplify(Complex(-27*sqrt(3)));
addressof(%); # The same as L[2]

 

You could rename all this way:
 

restart;
S1:=singular(1/sin(x),x);                    
S2:=singular(x/sin(x),x);
S3:=singular(x^2/sin(x),x);
indets(`union`(S1,S2,S3),`local`)=~_Z; 
assign(%);
`union`(S1,S2,S3); 
 # result  {x = Pi*_Z}

 

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