Maple Questions and Posts

These are Posts and Questions associated with the product, Maple

Hello. I want to solve a differential equation on a thickening grid. That is, to be solved first at N=10, then N=20, then N=40 and so on.
Tried the design
for N from 10 to 160 by 2*N do
but for Maple it is difficult. Do you have any ideas how to implement this type of cycle?

I do not understand why Maple sometimes shows singular solution to Clairaut ODE and sometimes not.

Clairaut ODE has the form y(x) = x y'(x) + G(x, y')

In the following ODE when I ask Maple to dsolve it as is, it does give singular solution. Next, when solving explicity for y(x) first, which will generate 2 ODE's, each is Clairaut ODE, then ask Maple to dsolve each, now Maple no longer gives the singular solution. But when I solve each one of these ODE's, I see that there is the singular solution there. It must be there, since this is Clairaut ODE and it has singular solution.

When I do PDEtools:-casesplit on each of the two ODE's generated by solving for y(x) first, I see the singular solution there.

The question is, why Maple dsolve does not show the singular solution in the second case? And how to make it show it? Or did I do something wrong?

restart;

Typesetting:-Settings(typesetprime=true):

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

x^2*(diff(y(x), x))^2-(1+2*x*y(x))*(diff(y(x), x))+1+y(x)^2 = 0

DEtools:-odeadvisor(ode)

[[_1st_order, _with_linear_symmetries], _rational, _Clairaut]

Vector([dsolve(ode,y(x))]); #now it shows singular solution (first one below)

Vector(3, {(1) = y(x) = (1/4)*(4*x^2-1)/x, (2) = y(x) = _C1*x-sqrt(_C1-1), (3) = y(x) = _C1*x+sqrt(_C1-1)})

PDEtools:-casesplit(ode)

`casesplit/ans`([(diff(y(x), x))^2 = (2*y(x)*(diff(y(x), x))*x+diff(y(x), x)-y(x)^2-1)/x^2], [2*(diff(y(x), x))*x^2-2*x*y(x)-1 <> 0]), `casesplit/ans`([y(x) = (1/4)*(4*x^2-1)/x], [])

ode:=convert(ode,D): #solve for y(x) first, this will generate 2 ODE's
sol:=[solve(ode,y(x))]:
odes:=Vector(map(z->y(x)=z,convert(sol,diff)))

Vector(2, {(1) = y(x) = (diff(y(x), x))*x+sqrt(diff(y(x), x)-1), (2) = y(x) = (diff(y(x), x))*x-sqrt(diff(y(x), x)-1)})

DEtools:-odeadvisor(odes[1]);
DEtools:-odeadvisor(odes[2]);

[[_1st_order, _with_linear_symmetries], _rational, _Clairaut]

[[_1st_order, _with_linear_symmetries], _rational, _Clairaut]

dsolve(odes[1],y(x)); #where is singular solution?

y(x) = _C1*x+(_C1-1)^(1/2)

dsolve(odes[2],y(x)); #where is singular solution?

y(x) = _C1*x-(_C1-1)^(1/2)

PDEtools:-casesplit(odes[1])

`casesplit/ans`([diff(y(x), x)-1 = (-(diff(y(x), x)-1)^(1/2)+y(x)-x)/x, diff(y(x), x) = (-(diff(y(x), x)-1)^(1/2)+y(x))/x], [1+2*x*(diff(y(x), x)-1)^(1/2) <> 0]), `casesplit/ans`([y(x) = (1/4)*(4*x^2-1)/x], [])

PDEtools:-casesplit(odes[2])

`casesplit/ans`([diff(y(x), x)-1 = ((diff(y(x), x)-1)^(1/2)+y(x)-x)/x, diff(y(x), x) = ((diff(y(x), x)-1)^(1/2)+y(x))/x], [2*x*(diff(y(x), x)-1)^(1/2)-1 <> 0]), `casesplit/ans`([y(x) = (1/4)*(4*x^2-1)/x], [])

 

 

Download missing_singular.mw

I would like to draw an ellipse by orthonal affinity of a circle. Thank you.

Dear Maple users

I wanted to export two long lists X and Y of data into two columns in Excel by using the Export command from the ExcelTools package. Data was however exported into rows! I converted the lists into vectors and I succeeded in doing it the way I wanted:

with(ExcelTools);
X1 := convert(X, Vector);
Y1 := convert(Y, Vector);
Export(X1, "data1.xlsx", 1, "A1");
Export(Y1, "data1.xlsx", 1, "B1");

But is it really necessary to convert to vectors in order to accomplish this task? 

Regards,

Erik


 

 

This is a problem that intrigues me.

I have been thinking about the distribution of n points, p[1], p[2], ..., p[n] in a scatterplot and I am interested in finding the location of a focal point, A, such that the distance from A to each neighboring point is a minimum. I have worked on a routine and this is given in the attached worksheet. In this example, n=20 and the position of A is determined to be [-3.25, 0.99].

Now - here's my question. Suppose I am interested in n-large and, instead of locating one focal point, A, I wish to obtain several,i.e. A, B, C, etc. such that the distances from each of these to their respective neighboring points is also a minimum.

Does any interested party know if this is possible to do and if so, can anyone suggest an approach or routine? If so, I'd be delighted to understand how to solve for this.

Thank you for reading!

Mapleprimes_Aug_17.mw

 

 

The commands

with(LinearAlgebra):interface(rtablesize=50);Matrix([[1,1],[4,5]]);

give the output


50
Matrix(2, 2, {(1, 1) = 1, (1, 2) = 1, (2, 1) = 4, (2, 2) = 5}, datatype = anything, storage = rectangular, order = Fortran_order, shape = [])

 

Why won't it display properly?

 

To help decide the type of ODE, I need a function which moves all derivatives to lhs and everything else to the rhs of the ODE. This way I can more easily analyze the ODE.

The ode will be only first order. Any term which contains  (y')^n, is to be moved to the lhs. Rest to rhs. Couple of examples will help illustrate the problem. This is done in code only, without looking at the screen. The dependent variable is always y and the independent variable is x.  So I need to turn the ODE to the form 

                    G(y',y'^2,y'^3,.....,y'^n)  = F(x,y)

I can find all derivatives in the ODE using

indets['flat'](ode,{`^`('identical'(diff(y(x),x)),'algebraic'),'identical'(diff(y(x),x))})

But not sure how this will help me do what I want. isolate does not help, since it only takes one term at a time. Collect also did not help for same reason.  If the ODE contains only ONE derivative, then it is easy to do. But the question is about how to do it for ODE which contains more than y' term of different powers as the examples below show.

What methods to do this in Maple? I looked at DEtools but so far did not see anything.

Example 1

 

restart;
ode:=3*diff(y(x),x)^2+diff(y(x),x)^3+sin(x)+y(x)=x*y(x)+x*diff(y(x),x);

3*(diff(y(x), x))^2+(diff(y(x), x))^3+sin(x)+y(x) = x*y(x)+x*(diff(y(x), x))

indets['flat'](ode,{`^`('identical'(diff(y(x),x)),'algebraic'),'identical'(diff(y(x),x))})

{(diff(y(x), x))^2, (diff(y(x), x))^3, diff(y(x), x)}

ode_wanted:= 3*diff(y(x),x)^2+diff(y(x),x)^3-x*diff(y(x),x)=-sin(x)-y(x)+x*y(x)

3*(diff(y(x), x))^2+(diff(y(x), x))^3-x*(diff(y(x), x)) = -sin(x)-y(x)+x*y(x)

Example 2

 

restart;
ode:=3*diff(y(x),x)^2+diff(y(x),x)=x*diff(y(x),x)+5;

3*(diff(y(x), x))^2+diff(y(x), x) = x*(diff(y(x), x))+5

indets['flat'](ode,{`^`('identical'(diff(y(x),x)),'algebraic'),'identical'(diff(y(x),x))})

{(diff(y(x), x))^2, diff(y(x), x)}

ode_wanted:= 3*diff(y(x),x)^2+diff(y(x),x)-x*diff(y(x),x)=5

3*(diff(y(x), x))^2+diff(y(x), x)-x*(diff(y(x), x)) = 5

 

 

Download lhs_rhs.mw

 

i want to solve the systems of diff equation what's the problem   0.mw

Dear Maple users

Some students have come to us to report, that something doesn't seem to work properly in Maple 2019.1 in Document Mode. And they seem to be right: writing an passive math formula by using Shift+F5 (the formula is gray, not blue), then using F5 to get out of that Math field and back into Text Mode. Using the Enter key to go to the next line: It doesn't work! The cursor stays in the same line. This behavior is new in Maple 2019. It worked properly in Maple 2018 and earlier. I assume it is not the intention? 

I know it can easily be dealt with by making a new Paragraph by using the shortcut Ctrl+Shift+J. I call the assumed bug 'severe' though, because it will severely delay the workflow for many students. They are used to deliver a document mixed with formulas (active or passive) and text. 

NB! I have tested it on several computers (Mac and Windows), and it doesn't work on any of them.

Regards,

Erik V.

Hi I have the following ln function that I want to differentiate wrt variable c:

 

I2 := -sqrt(-c^2+1)*ln(abs((.9*sqrt(-c^2+1)-c*sqrt(1-.9^2))/(-.9*sqrt(-c^2+1)-c*sqrt(1-.9^2))))

 

When I differentiate I obtain an expression that involves 

...abs(1, (.9*sqrt(-c^2+1)-.4358898944*c)/(-.9*sqrt(-c^2+1)-.4358898944*c))...

Why does it give a comma at the abs expression? how to get rid of that.

Seeking for fast approximate formulas to compute (a huge number of) quantiles of a Gaussian random variable (here the standard one, but its extension to any Gaussian RV is straightforward), I found a few of them in the Abramowitz and Stegun book, page 933, relations 26.2.22 and 26.2.23.
Each approximation model is expressed as a rational fraction, the second one being the more accurate.
Each model depends on (respectively 4 and 6) parameters that are estimated (I guess it was done this way) through a least-square-like method.

See here for an online access http://people.math.sfu.ca/~cbm/aands/page_933.htm.

These approximation, and specially the most accurate one (formula 26.2.23) seem to be still widely used today(1) (see for instance https://www.johndcook.com/blog/normal_cdf_inverse/ ).

As an amusement I decided to compute the best fit by using the Statistics:-NonLinearFit procedure and a sample of (probability, quantile) points where probability ranges in [0.5, 1-1/1000] (the range used in formulas 26.2.22 and 26.2.23 is (0, 0.5] but this is not a point).
Surprisingly Statistics:-NonLinearFit returned, for the two formulas, parameter estimations substantially different from the one given in the Abramowitz & Stegun's book. A reason could be that the points I used when I did the fits weren't the one they used (unfortunately they give no informations about this).

More interesting, whatever the formula I refitted,  NonLinearFit produced an approximation whose the absolute error was smaller by about two orders of magnitude to the onesprovided by Abramowitz and Stegun.
For instance they wrote that the most accurate formula (26.2.23) had an absolute approximation error less than 4.5*10-4 as I obtained a value around 10-6!

(1) To get an idea of the persistence of the use of the formula 26.2.23, just type the value 2.515517 of its parameter c[0] in any search engine.


In the plots below the gray rectangle refers to the region where the approximate ICDF is used for extrapolation.
 

restart:

with(Statistics):

cdf := unapply(evalf(CDF(Normal(0, 1), x)), x):
X   := [seq(0..5, 0.1)]:
A   := cdf~(X):
T  := alpha -> sqrt(-2*log(1-alpha)):
q  := Quantile~(Normal(0, 1), A):
Aq := convert([A,q], Matrix)^+:

r := 1:

J := z -> z - add(a__||k*z^k, k=0..r)/(1+add(b__||k*z^k, k=1..r+1)):


model  := J(T(alpha)):
NL_fit := unapply(NonlinearFit(model, Aq, alpha), alpha);


# these lines are for estimating the performances
B  := Sample(Uniform(0.5, 1), 10^4):
CodeTools:-Usage(Quantile~(Normal(0, 1), B)):
CodeTools:-Usage(Quantile~(Normal(0, 1), B, numeric)):
CodeTools:-Usage(NL_fit~(B)):
#-----------------------------------------------------
Y  := [seq(0..6, 0.01)]:
B  := cdf~(Y):
R1 := Quantile~(Normal(0, 1), B, numeric):
R2 := NL_fit~(B):

plots:-display(
  ScatterPlot(R1, log[10]~(abs~(R2-~R1)), legend=mmcdara, color=red, gridlines=true, size=[700, 400]),
  plottools:-rectangle([max(X), log[10]~(min(abs~(R2-~R1)))], [max(Y), log[10]~(max(abs~(R2-~R1)))], color=gray, transparency=0.6)
);

proc (alpha) options operator, arrow; (-2*ln(1-alpha))^(1/2)-(HFloat(2.5454311687345044)+HFloat(0.8058592540791468)*(-2*ln(1-alpha))^(1/2))/(1+HFloat(1.4689746699940707)*(-2*ln(1-alpha))^(1/2)-HFloat(0.34455942407858625)*ln(1-alpha)) end proc

 

memory used=170.31MiB, alloc change=76.01MiB, cpu time=3.06s, real time=3.05s, gc time=54.87ms

memory used=171.59MiB, alloc change=256.00MiB, cpu time=3.12s, real time=3.03s, gc time=154.77ms

memory used=8.24MiB, alloc change=0 bytes, cpu time=95.00ms, real time=95.00ms, gc time=0ns

 

 

r := 2:

 
J := z -> z - add(a__||k*z^k, k=0..r)/(1+add(b__||k*z^k, k=1..r+1)):


model  := J(T(alpha)):
NL_fit := unapply(NonlinearFit(model, Aq, alpha), alpha);


# these lines are for estimating the performances
B  := Sample(Uniform(0.5, 1), 10^4):
CodeTools:-Usage(Quantile~(Normal(0, 1), B)):
CodeTools:-Usage(Quantile~(Normal(0, 1), B, numeric)):
CodeTools:-Usage(NL_fit~(B)):
#-----------------------------------------------------


Y  := [seq(0..6, 0.01)]:
B  := cdf~(Y):
R1 := Quantile~(Normal(0, 1), B, numeric):
R2 := NL_fit~(B):

plots:-display(
  ScatterPlot(R1, log[10]~(abs~(R2-~R1)), legend=mmcdara, color=red, gridlines=true, size=[700, 400]),
  plottools:-rectangle([max(X), log[10]~(min(abs~(R2-~R1)))], [max(Y), log[10]~(max(abs~(R2-~R1)))], color=gray, transparency=0.6)
);

proc (alpha) options operator, arrow; (-2*ln(1-alpha))^(1/2)-(HFloat(2.9637294443959394)+HFloat(4.527738737327481)*(-2*ln(1-alpha))^(1/2)-HFloat(0.9571637188191973)*ln(1-alpha))/(1+HFloat(3.472400103322335)*(-2*ln(1-alpha))^(1/2)-HFloat(3.426536241250657)*ln(1-alpha)+HFloat(0.08875278252087411)*(-2*ln(1-alpha))^(3/2)) end proc

 

memory used=170.09MiB, alloc change=32.00MiB, cpu time=3.29s, real time=3.11s, gc time=268.60ms

memory used=170.85MiB, alloc change=0 bytes, cpu time=3.23s, real time=3.10s, gc time=201.52ms
memory used=10.76MiB, alloc change=0 bytes, cpu time=127.00ms, real time=127.00ms, gc time=0ns

 

 

# Optimized "r=2" computation

z_fit := simplify(subs(alpha=-exp(-(1/2)*z^2)+1, NL_fit(alpha))) assuming z > 0:
z_fit := unapply(convert~(%, horner), z);

p := proc(alpha)
  local z:
  z := sqrt(-2*log(1-alpha)):
  z_fit(z):
end proc:

R3 := CodeTools:-Usage(p~(B)):

plots:-display(
  ScatterPlot(R1, log[10]~(abs~(R2-~R1)), legend=mmcdara, color=red, gridlines=true, size=[700, 400]),
  plottools:-rectangle([max(X), log[10]~(min(abs~(R2-~R1)))], [max(Y), log[10]~(max(abs~(R2-~R1)))], color=gray, transparency=0.6)
);

proc (z) options operator, arrow; (-2.963729444+(-3.527738737+(2.993818244+(1.713268121+0.8875278252e-1*z)*z)*z)*z)/(1.+(3.472400103+(1.713268121+0.8875278252e-1*z)*z)*z) end proc

 

memory used=1.67MiB, alloc change=0 bytes, cpu time=14.00ms, real time=15.00ms, gc time=0ns

 

 


AS stands for Abramowith & Stegun

J_AS := unapply(normal(eval(J(t), [a__0=2.515517, a__1=0.802853, a__2=0.010328, b__1=1.432788, b__2=0.189269, b__3=0.001308])), t):
J_AS(t);


# for comparison:

print():
z_fit := simplify(subs(alpha=-exp(-(1/2)*z^2)+1, NL_fit(alpha))) assuming z > 0:
map(sort, %, z);

plot([z_fit(z), J_AS(z)], z=0.5..1, color=[blue, red], legend=[mmcdara, Abramowitz_Stegun], gridlines=true);

print():
R2_AS := CodeTools:-Usage(J_AS~(T~(B))):
print():


plots:-display(
  ScatterPlot(R1, log[10]~(abs~(R2_AS-~R1)), legend=Abramowitz_Stegun, gridlines=true, size=[700, 400]),
  ScatterPlot(R1, log[10]~(abs~(R2-~R1)), legend=mmcdara, color=red),
  plottools:-rectangle([max(X), log[10]~(min(abs~(R2-~R1)))], [max(Y), log[10]~(max(abs~(R2-~R1)))], color=gray, transparency=0.6)
);

(0.1308000000e-2*t^4+.1892690000*t^3+1.422460000*t^2+.1971470000*t-2.515517000)/(0.1308000000e-2*t^3+.1892690000*t^2+1.432788000*t+1.)

 

 

(0.8875278252e-1*z^4+1.713268121*z^3+2.993818244*z^2-3.527738737*z-2.963729444)/(0.8875278252e-1*z^3+1.713268121*z^2+3.472400103*z+1.)

 

 

 

memory used=2.92MiB, alloc change=0 bytes, cpu time=25.00ms, real time=25.00ms, gc time=0ns

 

 

 


 

Download InverseNormalCDF.mw

 

 

Although the graph of a parametrized surface can be viewed and manipulated on the computer screen as a surface in 3D, it is not quite suitable for printing on a 3D printer since such a surface has zero thickness, and thus it does not correspond to physical object.

To produce a 3D printout of a surface, it needs to be endowed with some "thickness".  To do that, we move every point from the surface in the direction of that point's nomral vector by the amount ±T/2, where T is the desired thickness.  The locus of the points thus obtained forms a thin shell of thickness T around the original surface, thus making it into a proper solid. The result then may be saved into a file in the STL format and be sent to a 3D printner for reproduction.

The worksheet attached to this post provides a facility for translating a parametrized surface into an STL file.  It also provides a command for viewing the thickened object on the screen.  The details are documented within that worksheet.

Here are a few samples.  Each sample is shown twice—one as it appears within Maple, and another as viewed by loading the STL file into MeshLab which is a free mesh viewing/manipulation software.

 

Here is the worksheet that produced these:  thicken.mw

 

 

i have a question about sets. how can i keep set members in order of addition not the defualt maple ordering.
in maple help, i saw the command setsort=0..8 but i do not know how to use it.
consider exmaple below:


 

restart

L:={}:

for j in [3,5,6,1] do
birth:=j:
L:=`union`(L ,{j});
od;

3

 

{3}

 

5

 

{3, 5}

 

6

 

{3, 5, 6}

 

1

 

{1, 3, 5, 6}

(1)

 


how can i keep L as order of addition? L={3,5,6,1}. thanks in advance

Download setoder.mw

Hi I am using Maple 18 and am trying to make my looping work backward.

 

with(Student[Calculus1]);
with(linalg);
with(orthopoly);
Digits := 20;
n, K := 5, 4;
for k from 0 to n do a[k] := (int((T(2, t)-T(0, t))*T(k, t)/sqrt(-t^2+1), t = -1 .. 1))/Pi end do

First I use the above command to generate a0, a1, a2, a3, a4 and a5. Now i want to do a back substitute using the following recurrence relation to determine the values of b5, b4, b3, b2, b1 and b0.

b[n] := 0; (to find b5)

b[n-1] := a[n]; (to find b4)

(to find b3, b2, b1, b0)

Can someone assist me with the backward looping please. I does not seem to work.

 

 

I calculate volume of tetrahedron SABC by this way. Now, I want to find one option  of coordinates of all vertices of this tetrahedron. How can I get the coordinates of the points S, A, B, C so that coordinates  center of circumcircle of the triangle ABC is (0,0,0)?

 

restart:
SA := 3: 
SB := 5: 
SC := 7:
AB := 3:
BC := 4: 
AC := 5:
V := (1/12)*sqrt(SA^2*BC^2*(AB^2+AC^2-BC^2-SA^2+SB^2+SC^2)+SB^2*AC^2*(AB^2-AC^2+BC^2+SA^2-SB^2+SC^2)+SC^2*AB^2*(-AB^2+AC^2+BC^2+SA^2+SB^2-SC^2)-SA^2*AB^2*SB^2-SB^2*BC^2*SC^2-SA^2*AC^2*SC^2-AB^2*BC^2*AC^2);

From here https://www.mapleprimes.com/questions/227573-How-To-Get-Coordinates-Of-Triangle-ABC
I find coordinates of the point A, B, C. Now I tried
 

with(geom3d):
point(A, -3/2, -2, 0):
point(B, 3/2, -2, 0):
point(C, -3/2, 2, 0):
point(S, x, y, z):
solve([distance(S, A) = 3, distance(S, B) = 5, distance(S, C) = 7], [x, y, z], explicit);
with(geom3d):
point(A, -3/2, -2, 0):
point(B, 3/2, -2, 0):
point(C, -3/2, 2, 0):
point(S, x, y, z):
solve([distance(S, A) = 5, distance(S, B) = 3, distance(S, C) = 7], [x, y, z], explicit);

Use the above code, I obtain the result.

1 2 3 4 5 6 7 Last Page 2 of 1574