## 1350 Reputation

15 years, 11 days
University of Twente (retired)
Enschede, Netherlands

My "website" consists of a Maple Manual in Dutch

## Keeping the curve...

You will see that the circular curve disappears in the last frame. This has to do with rounding errors. The last sqrt is a complex number with a small imaginary part. You can avoid this by replacing the third command by

`C := animatecurve( [1/2*p,1/2*Re(sqrt(4-p^2)), p = 0..2], frames=50 ):`

## Multiple solutions...

@GrecoRoman Firstly, in Error.mw the formal parameter is data, not x. Secondly, the solve command in EuclideanRegression has multiple solutions as can be seen by adding a print command:

Moreover, in the last unapply statement the y must be replaced by x

In Euclidean.mw the third solution is apparantly the solution you want.

## Multiple solutions...

@GrecoRoman Firstly, in Error.mw the formal parameter is data, not x. Secondly, the solve command in EuclideanRegression has multiple solutions as can be seen by adding a print command:

Moreover, in the last unapply statement the y must be replaced by x

In Euclidean.mw the third solution is apparantly the solution you want.

## Right!...

@kjetil87  , this is the right approach.
The sum command tries to calculate a symbolic sum if the range contains parameters or infinity. If the range contains only numbers you can better use the add command. When calculating a symbolic sum, only the generic values of the variables are taken into account; eventual singularities you have to find out for yourself. In your example it is easily seen that only p1=0 and p1=0 lead to special cases that have to do with "00 = 1" So the best thing you can do is indeed

`restart;A := p1^abs(l)*p2^abs(l-2)*exp(-I*2*Pi*f*l*T):sum( A, l=1..infinity); # generic case                               p1                               ----------------------------                  p2 (exp(2 I Pi f T) - p1 p2)S := add(A, l=1..2) + sum(A, l=3..infinity): # special caseseval(S,p2=0);                        2                                       p1  exp(-4 I Pi f T)eval(S,p1=0);                               0`

## Right!...

@kjetil87  , this is the right approach.
The sum command tries to calculate a symbolic sum if the range contains parameters or infinity. If the range contains only numbers you can better use the add command. When calculating a symbolic sum, only the generic values of the variables are taken into account; eventual singularities you have to find out for yourself. In your example it is easily seen that only p1=0 and p1=0 lead to special cases that have to do with "00 = 1" So the best thing you can do is indeed

`restart;A := p1^abs(l)*p2^abs(l-2)*exp(-I*2*Pi*f*l*T):sum( A, l=1..infinity); # generic case                               p1                               ----------------------------                  p2 (exp(2 I Pi f T) - p1 p2)S := add(A, l=1..2) + sum(A, l=3..infinity): # special caseseval(S,p2=0);                        2                                       p1  exp(-4 I Pi f T)eval(S,p1=0);                               0`

## evalb...

@jam to test if the left hand side of an equation equals the right hand side you can use ?evalb (after simplify).

## evalb...

@jam to test if the left hand side of an equation equals the right hand side you can use ?evalb (after simplify).

The function that you want to plot is only defined on the circular disk x2 + y2 ≤ 25.

Therefore you can use polar coordinates for x and y to obtain constant ranges for the parameters:

`plot3d( subs( {x=r*cos(t),y=r*sin(t)}, [x,y,sqrt(25-x^2-y^2)] ), r=0..5, t=0..2*Pi, axes=boxed );`

Another possibility is

` plot3d( sqrt(25-x^2-y^2), x=-5..5, y=-sqrt(25-x^2)..sqrt(25-x^2), axes=boxed );`

But this also misses a few points.

The function that you want to plot is only defined on the circular disk x2 + y2 ≤ 25.

Therefore you can use polar coordinates for x and y to obtain constant ranges for the parameters:

`plot3d( subs( {x=r*cos(t),y=r*sin(t)}, [x,y,sqrt(25-x^2-y^2)] ), r=0..5, t=0..2*Pi, axes=boxed );`

Another possibility is

` plot3d( sqrt(25-x^2-y^2), x=-5..5, y=-sqrt(25-x^2)..sqrt(25-x^2), axes=boxed );`

But this also misses a few points.

## Indexing function...

@acer , I also tried to find an indexing function, but I didn't succeed.
But your quick solution doesn't work:

`restart;m := LinearAlgebra:-RandomMatrix(4,6):M := ArrayTools:-UpperTriangle(m);V := Vector( convert( select( x -> x<>0, M ), list ) ):G:=(m,n,W) -> Matrix(m,n,(i,j)->`if`(i>j,0,W[(i-1)*n-((i-1)*(i)/2)+j])):G(4,6,V);`

## Indexing function...

@acer , I also tried to find an indexing function, but I didn't succeed.
But your quick solution doesn't work:

`restart;m := LinearAlgebra:-RandomMatrix(4,6):M := ArrayTools:-UpperTriangle(m);V := Vector( convert( select( x -> x<>0, M ), list ) ):G:=(m,n,W) -> Matrix(m,n,(i,j)->`if`(i>j,0,W[(i-1)*n-((i-1)*(i)/2)+j])):G(4,6,V);`

## Try this...

Casper, I'm not sure if I have all the indices right. But I think this may help you further.

`restart:para := proc(C,new)   global p;    # new = 0 or 1    seq(log(p[j,C]/(1-p[j,C]))=mu+tau[j]+eta[C] + new*tau[j]*eta[C],j=2..4)end proc:getpVec:=proc(k,c,new)  local i,j, parC;  global p, wp;  for i to c do    solve({para(i,new)},[seq(p[j,i],j=2..(k+1))]);    parC[i] := Vector(k, [rhs~(op(%))] );  end do;  add(wp[i]*parC[i],i=1..(c-1))+(1-add(wp[i],i=1..(c-1)))*parC[c];end proc:getpVec(3,3,1);`

## Try this...

Casper, I'm not sure if I have all the indices right. But I think this may help you further.

`restart:para := proc(C,new)   global p;    # new = 0 or 1    seq(log(p[j,C]/(1-p[j,C]))=mu+tau[j]+eta[C] + new*tau[j]*eta[C],j=2..4)end proc:getpVec:=proc(k,c,new)  local i,j, parC;  global p, wp;  for i to c do    solve({para(i,new)},[seq(p[j,i],j=2..(k+1))]);    parC[i] := Vector(k, [rhs~(op(%))] );  end do;  add(wp[i]*parC[i],i=1..(c-1))+(1-add(wp[i],i=1..(c-1)))*parC[c];end proc:getpVec(3,3,1);`

## Reconstruction of M...

...so I make a procedure to do this:

`RestoreMatrix := proc(r,c,V)local M, i,j,k;M := Matrix(r,c):k := 1:for i to r do   for j to i do     M[j,i] := V[k]:     k := k+1   end doend do;for i from r+1 to c do   for j to r do     M[j,i] := V[k];     k := k+1   end doend do;return Mend proc;RestoreMatrix(4,6,V);`

## Reconstruction of M...

...so I make a procedure to do this:

`RestoreMatrix := proc(r,c,V)local M, i,j,k;M := Matrix(r,c):k := 1:for i to r do   for j to i do     M[j,i] := V[k]:     k := k+1   end doend do;for i from r+1 to c do   for j to r do     M[j,i] := V[k];     k := k+1   end doend do;return Mend proc;RestoreMatrix(4,6,V);`
 5 6 7 8 9 10 11 Page 7 of 11
﻿