## 26 Badges

14 years, 206 days

## CurveFitting:-LeastSquares...

Because these 3 points are located in the same plane (the plane  yOz ), then you can use the  CurveFitting:-LeastSquares  command. We use a 1st degree polynomial  z=a*y+b . Increasing the degree gives nothing, because you have only 3 points, and two of them lie on the same vertical. A polynomial of any degree  z=P(y)  found by the least squares method will pass through point  C  and the middle of the segment  AB :

 > restart; with(plots): A:=[0,0,0]: B:=[0,0,1]: C:=[0,1,0]: eq:=CurveFitting:-LeastSquares([A[2..3],B[2..3],C[2..3]], y); P1:=spacecurve([0,y,eq], y=-1..2, color=green, thickness=2): P2:=pointplot3d([A,B,C],color=[blue,red,black], symbol=cross, symbolsize=40, thickness=2): P3:=textplot3d([[A[],"A"],[B[],"B"],[C[],"C"]], font=[TIMES,14], align={above,left}): display(P1,P2,P3, labels=[x,y,z], view=[(-1..1)\$3], orientation=[125,65]);
 >

Download 3_points.mw

## Adjustment and plotting...

1. I made some corrections to your worksheet that Thomas Richard mentioned already.

2. I am not quite sure what exactly you want to plot. For each value of  y , equations   and  b  define some lines on the plane. Below on the first graph these lines are plotted for  y = 1 . On the second graph, using the  Explore  command, you can trace how these graphs change with the change in the variable  y .

 > restart; # z=0.001, 0.002  and  r=0.1-0.3 para := z=0.001, r=0.1: eq := 1-((z/x^2)+(1-r)/(x-y)^2+r/(x+((1-r)/r)*y)^2)=0: eq1 := subs({x=p+q*I,para},eq); a := evalc(Re(eq1)); b := evalc(Im(eq1)); plots:-implicitplot(eval([a,b],y=1), p=-0.3..2.3, q=-0.5..0.5, color=[red,blue], thickness=2, gridrefine=4, scaling=constrained, size=[1000,300], axes=box, legend=["a","b"], title = "The plots of a and b for y=1", titlefont=[TIMES,16]); Explore(plots:-implicitplot([a,b], p=-5..3, q=-0.5..0.5, color=[red,blue], thickness=2, gridrefine=4, scaling=constrained, size=[800,300], axes=box, legend=["a","b"], title = "The plots of a and b for y=0..2", titlefont=[TIMES,16]), y=0...2.);
 >

Download eq_(1)_new.mw

## plots:-contourplot...

Should be the  plots:-contourplot  command instead of the  contourplot one,  because it is from a pluggable package plots.

## plot3d...

Remove the extra letter from the command name. Should be  plot3d   instead of  plotd3d .

## Solution...

 > restart;
 > f:=(x,y)->x^4-3*x^2-2*y^3+3*y+1/2*x*y;
 (1)

Find the critical points:

 > CriticalPoints:= [solve( {diff(f(x,y),x)=0,diff(f(x,y),y)=0}, {x,y}, explicit)]: evalf(%);
 (2)

Find the second partial derivatives:

 > A:=diff(f(x,y),x,x); B:=diff(f(x,y),x,y); C:=diff(f(x,y),y,y);
 (3)

In a loop, using the quadratic form  A*C-B^2, examine each of the 6 critical points:

 > for k from 1 to 6 do S:=evalf(CriticalPoints[k]); Q:=eval(A*C-B^2,S); A0:=eval(A,S); f0:=eval(f(x,y),S); if Q>0 then if A0<0 then R[k]:=[Maximum,S,f0] else R[k]:=[Minimum,S,f0] fi;fi; if Q<0 then R[k]:=[Saddle,S,f0] fi; od: Matrix([[Type_of_point,Point,Value_of_function],seq(R[k],k=1..6)]);
 (4)

Plotting:

 > P1:=plot3d(f(x,y),x=-2..2,y=-2..2, view=-6..6): P2:=plots:-pointplot3d([seq(eval([x,y,z],[R[k][2][],z=R[k][3]]),k=1..6)],color=red, symbol=solidsphere, symbolsize=18): plots:-display(P1,P2, orientation=[60,60]);
 >

The plot clearly shows: one maximum, two minimums and three saddles

Download zadelpunten_procedure_new.mw

## seq...

Use the  seq  command for this.
Example:

```restart;
L:=[0, Pi/6, Pi/4, Pi/3, Pi/2];
[seq([p,sin(p)], p=L)];```

## PolynomialTools:-Split, convert(..., ra...

```restart;
P:=x^4 - x^2 + x - 1:
A:=PolynomialTools:-Split(P, x);
convert(A, radical);
simplify(expand(%));  # A check
```

## Solution...

1. For  isolve a system should be specified as a set (not a list).
2. You need to specify exactly what you mean by the smallest positive solution, because the solution will not be one number, but an ordered set of numbers (a list). For example, which is less: [2, 4] or [1, 5] ?

In the solution below, the solutions are sorted by the minimum sum of the squares of all coordinates.

```restart;
sys:= {a[1] = 2 * a[5], a[1] = a[4], 4 * a[1] + 2 * a[3] = a[6] + 2 * a[7], 2 * a[2] + 2 * a[3] = 2 * a[6], a[2] = a[4] + a[5], 2 * a[2] -a[1] = 2 * a[4]};
indets(sys);
n:=nops(%);
isolve(sys);
assign(%);
solve({seq(a[i]>0, i=1..n)});
[seq(seq([seq(a[i],i=1..n)],_Z2=floor(5*_Z1/2)+1..10),_Z1=1..10)];
sort(%, key = (x->add(x[i]^2,i=1..n)));
```

## Round-off error...

I think this contradiction is just a rounding error. Maple uses 10 digits by default. Use precise arithmetic (or make  Digits  bigger):

```convert(-24.*y^6*(y-1.)^7*(y+1.)^7*(4.+6.*y^8-22.*y^6+30.*y^4-18.*y^2)/(4.*y^18-28.*y^16+84.*y^14+4.-136.*y^12+116.*y^10-24.*y^8-52.*y^6+56.*y^4-24.*y^2)^2, fraction);
subs(y=1+1/100000, %);
evalf(%);
```

## Solution...

First, we find the equation of the normal at an arbitrary point  P(x1,x1^2)  on the parabola. This normal intersects the axis  Oy  at a point  Q(0,y0) . Imposing the condition that the distance  PQ  = 1, we first find x1, and then the point  Q (the center of desired circle).

```restart;
f:=x->x^2:
g:=x->f(x1)-1/D(f)(x1)*(x-x1):
y0:=g(0):
P,Q:=[x1,f(x1)],[0,y0];
solve((P[1]-Q[1])^2+(P[2]-Q[2])^2=1);
x1:=%[1]:
Q;

# Visualization
A:=plot([f(x),[t,g(t),t=Q[1]..P[1]]], x=-2..2):
B:=plot([P,Q],style=point):
C:=plots:-textplot([[P[],"P"],[Q[],"Q"]], align=right, font=[TIMES,14]):
c:=plottools:-circle(Q,1):
plots:-display(A,B,C,c, scaling=constrained);
```

## A way...

Yesterday I answered the same question only asked in another post (this post disappeared for some reason). Here's my answer:

```restart;
A__left:=1/C__AO-1/C__AA;
f:=1/C__AO;
A__right:=``(f)*(expand(A__left/f));
expand(A__right); # Check
```

## A way...

```restart;
A__left:=1/C__AO-1/C__AA;
f:=1/C__AO;
A__right:=``(f)*(expand(A__left/f));
expand(A__right); # Check
```

## evalc...

@jrive  For arbitrary real parameters the last line should be

Pabang := evalc(argument(Pab));

Returns

From help  "For real arguments x, y, the two-argument function arctan(y, x), computes the principal value of the argument of the complex number
x + I y
, so
-Pi < arctan(y, x) and arctan(y, x) <= Pi"

## Implementation with Maple...

The code below implements this process in Maple. As already pointed out by vv, on each new launch, the numbers  N  may be different.

```restart;
randomize():
n:=2012:
S:={seq([i,i], i=1..n)}:
for step from 1 to n-2 do
t:=combinat:-randcomb(S, 2);
S:=(S minus t) union {[n+step,2*(t[1][2]+t[2][2])]};
od:
N:=S[1][2]+S[2][2];
N/1000;
```

## The link...

 First 11 12 13 14 15 16 17 Last Page 13 of 274
﻿