## 13613 Reputation

19 years, 229 days

## Statistics:-Fit...

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

## Move the initial R from 0 and right...

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

## You need different names...

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 and frontend...

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],[{`+`,`*`,`=`}]);

## subsindets...

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

## Use a name or a string...

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.

## A better way...

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.

## What do you think is wrong?...

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)]
interface(verboseproc=3);

Continuing with L:

L := [sqrt(3)*I, -27*I*sqrt(3)]
LL:=convert~(L,list);
###
###
###

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

Continuing a little further:

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

## Renaming all _Zn...

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}

## Try style=point...

Try this:

plot([cos(t), sin(t), t = 0 .. 2*Pi*10000],adaptive=false);
plot([cos(t), sin(t), t = 0 .. 2*Pi*10000],adaptive=false,style=point);
plot([cos(t), sin(t), t = 0 .. 2*Pi*10000],style=point);

The result from the first connects the points (those seen when using style=point, and adaptive = false).
The third generates many more points.
Try increasing the number of points in the first:

plot([cos(t), sin(t), t = 0 .. 2*Pi*10000],adaptive=false,numpoints=500);

With numpoints=500 you get close to seeing what is going on.

PS. So how many points are chosen by the last plot?
You can find out by doing:

A:=plot([cos(t), sin(t), t = 0 .. 2*Pi*10000]);
op(1,A); #5445 points (only)

To get the same data you can also do:

plottools:-getdata(A);

PPS.  I'm very fond of odeplot. Here is a plot made by odeplot:

restart;
ode:=diff(x(t),t,t)+x(t)=0;
sol:=dsolve({ode,x(0)=1,D(x)(0)=0}); # first argument to you plot
-diff(sol,t); # second argument to your plot
res:=dsolve({ode,x(0)=1,D(x)(0)=0},numeric,maxfun=0); # Numerical solution
plots:-odeplot(res,[x(t),-diff(x(t),t)],0..2*Pi);
plots:-odeplot(res,[x(t),-diff(x(t),t)],0..2*Pi*10000,numpoints=10000);

Another version that doesn't involve using an ode particular yo the plot you want to make:

restart;
sys:={diff(x(t),t)=0,x(0)=0,diff(y(t),t)=0,y(0)=0};
dsolve(sys); # As intended
resS:=dsolve(sys,numeric,maxfun=10^6);
plots:-odeplot(resS,[x(t)+cos(t),y(t)+sin(t)],0..2*Pi,labels=["",""]);
plots:-odeplot(resS,[x(t)+cos(t),y(t)+sin(t)],0..2*Pi*10000,numpoints=10000,labels=["",""]);

The problem is the expansion that occurs.
Here is an ad hoc solution:

restart;
integrand:=(2*x^2022+1)/(x^2023+x);

num:=numer(integrand);
den:=denom(integrand);
igt:=(num-x^2022)/den + x^2022/den;
map(int,igt,x);

integrand2:=(2*x^n+1)/(x^m+x);
int(integrand2,x);
eval(%,{n=2022,m=2023});

Generalizing further:

integrand3:=(a*x^n+b)/(c*x^m+d*x);
int(integrand3,x) assuming a>0,b>0,c>0,d>0;
eval(%,{a=2,b=1,c=1,d=1,n=2022,m=2023});

## Another way...

Instead of the version given by nm, you could do:

f := unapply(piecewise(0 < x and x < a, k1, a < x and x < b, k2), a,x);

simplify(f(a,x)) assuming 0 < a, 0 < x and x < a;
simplify(f(b,y)) assuming 0 < b, 0 < y and y < b;nn

Both return k1.
The difference between the two versions of simplify with assuming/assume, is that in your version any assumptions are placed on the expression before any evaluation takes place. But then your function fa(x) doesn't see a inside the procedure fa, and so no assumption is therefore placed on a.
In the version given by nm fa(x) is evaluated first and assumptions are made.
In my version above the assumptions are placed on the arguments given right away, in my two examples a,x and b,y.

## EQ3 is the problem...

Change EQ3 to:

EQ3 := c[11] + c[12] + b^2*Q[0, 1]*sin(n*Pi*y[0, 1]/b)*(exp(-n*Pi*x[0, 1]/b) - exp(n*Pi*x[0, 1]/b))/(n*Pi) = lambda* (c[21] + c[22]);

It was:

EQ3 := c[11] + c[12] + b^2*Q[0, 1]*sin(n*Pi*y[0, 1]/b)*(exp(-n*Pi*x[0, 1]/b) - exp(n*Pi*x[0, 1]/b))/(n*Pi) = lambda . (c[21] + c[22]);
# Notice the right hand side. The dot should be *

## Remove eval...

If you remove eval from the definition of h then you get:

f();
g says: z evaluates to 2.
h says: z evaluates to z.

You set z:=2;  in f and write h('z');

Inside the procedure f this evaluates to h(z) only, because of the unevaluation quotes.

That doesn't mean that h is being handed a variable without any value. It just need evaluation, which it receives in your version of h.

Consider this simple example:

restart;
p:=proc() local z; z:=2; 'z' end proc;
p(); # z
%; # Now full evaluation to 2

## applyrule is limited...

I think that applyrule isn't meant to be all that ambitious. A quote from the help page for applyrule:

"It is more powerful than the command subs, but does not do mathematical transformations as algsubs does."

The only thing I could come up with is this:

restart;
rule1:=A::`&+`(anything,anything)=bbbb;
applyrule(rule1,[a+b,c+d,a+b+c]);
rule2:=A::`&+`(anything,anything)=''`*`(op(A))'';  # Notice double unevaluation quotes
applyrule(rule2,[a+b,c+d,a+b+c]);

Note: If you change the type in 'the_rule" to symbol, you get a+b:

restart;

the_rule:=A::symbol+B::symbol=A*B;

applyrule(the_rule,a+b); # result a+b, NOT a*b

Now it gets weird because if you change to 'name' then you get what you want a*b:

restart;

the_rule:=A::name+B::name=A*B;

applyrule(the_rule,a+b); # Result a*b

So my conclusion is that applyrule is rather limited and when used outside of its intended use, it gives you strange results. Maybe an example of garbage in, garbage out.
The procedure has a copyright from 2005.

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