mmcdara

7344 Reputation

18 Badges

8 years, 289 days

MaplePrimes Activity


These are answers submitted by mmcdara

Is this what you want?
(assumming derivatives with respect to x and y do commute)
 

restart:

alias(V=V(x, y)):
alias(alpha=alpha(x,y)):
alias(beta=beta(x,y)):

cond_1 := diff(V, x)=V*alpha:
cond_2 := diff(V, y)=V*beta:

eval(diff(V, x$2), cond_1);
eval(diff(V, x$2), cond_2);

(diff(V, x))*alpha+V*(diff(alpha, x))

 

diff(diff(V, x), x)

(1)

U := copy(V):
for k from 1 to 2 do
  U := eval(diff(U, y), cond_2);
end do;

V*beta

 

V*beta^2+V*(diff(beta, y))

(2)

der := proc(r, s)
  local fx, fy:

  fx := proc(W, r)
    local U, k:
    U := copy(W);
    for k from 1 to r do
      U := eval(diff(U, y), cond_2);
    end do;
    return U:
  end proc:

  fy := proc(W, s)
    local U, k:
    U := copy(W);
    for k from 1 to s do
      U := eval(diff(U, y), cond_2);
    end do;
    return U:
  end proc:

  if r=0 then
    if s=0 then
      V
    else
      Diff(V, y$s) = fy(V, s)
    end if:
  else
    if s=0 then
      Diff(V, y$r) = fx(V, r)
    else
      Diff(Diff(V, x$r), y$s) = fx(fy(V, s), r); # or  fy(fx(V, r), s) assuming derivatives commute
    end if:
  end if:
end proc

proc (r, s) local fx, fy; fx := proc (W, r) local U, k; U := copy(W); for k to r do U := eval(diff(U, y), cond_2) end do; return U end proc; fy := proc (W, s) local U, k; U := copy(W); for k to s do U := eval(diff(U, y), cond_2) end do; return U end proc; if r = 0 then if s = 0 then V else Diff(V, `$`(y, s)) = fy(V, s) end if else if s = 0 then Diff(V, `$`(y, r)) = fx(V, r) else Diff(Diff(V, `$`(x, r)), `$`(y, s)) = fx(fy(V, s), r) end if end if end proc

(3)

printlevel := 2:
for r from 0 to 2 do
  for s from 0 to 2 do
    der(r, s)
  end do:
end do;

V

 

Diff(V, y) = V*beta

 

Diff(V, y, y) = V*beta^2+V*(diff(beta, y))

 

Diff(V, y) = V*beta

 

Diff(Diff(V, x), y) = V*beta^2+V*(diff(beta, y))

 

Diff(Diff(V, x), y, y) = V*beta^3+3*V*beta*(diff(beta, y))+V*(diff(diff(beta, y), y))

 

Diff(V, y, y) = V*beta^2+V*(diff(beta, y))

 

Diff(Diff(V, x, x), y) = V*beta^3+3*V*beta*(diff(beta, y))+V*(diff(diff(beta, y), y))

 

Diff(Diff(V, x, x), y, y) = V*beta^4+6*V*beta^2*(diff(beta, y))+3*V*(diff(beta, y))^2+4*V*beta*(diff(diff(beta, y), y))+V*(diff(diff(diff(beta, y), y), y))

(4)

 


 

Download derivatives.mw

 

Your function is extremely simple for it can be rewritten this way

ff = C + f (xx[1]) + ... f (xx[N]) 

where C is a constant and f(x) = x^2 - 10 * cos(2*Pi*x)  (A and B are two constants)
Note that the hessien of ff iis thus a diagonal matrix

Each component of the gradient of ff is defined by G := u -> 2*u + 20*Pi * sin(2*Pi*u) .
Thus each component of the gradient of ff  is null for exactly the same values of xx[n].
In the range xx[n]=-5..5 (for each n) the function u -> 2*u + 20*Pi * sin(2*Pi*u)  has 21 zeroes.
To find their values use this pice of code
Z := NULL:
z := -5:
k := 1;
while k <= 21 do
  z := RootFinding:-NextZero(u->2*u+20*Pi*sin(2*Pi*u), z);
  Z := Z, z:
  k := k+1:
end do:
Z;


The 2nd derivative of u -> 2*u - 10 * cos(2*Pi*u) is  the function H := u -> 2 + 20*Pi * sin(2*Pi*u)
Note that the hessien of ff is a diagonal matrix whose diagonal elements are just H(xx[n])
The values of H(u) for each valus of Z are calculated by
Hval := map(u -> 2+20*Pi*sin(2*Pi*u), [Z])

Among these 21 values, 13 of them are poitive :
numelems(select[flatten](`>`, Hval, 0))

Then ff has a minimum if all the components (1000) of its hessian are negative at the points where the gradient of ff is null, and has a maximum if all the components of its hessian are positve at the points where the gradient of ff is null,
As 13 values of Hval are positive, one might expect 13^1000 minima (and 8^1000 maxima) for ff.
The number of saddle points is 21^1000 - 13^1000 - 8^1000

The largest value of ff at a point where its gradient is 0, is C + 1000 * V where V is the highest value of 2*u - 10 * cos(2*Pi*u) (u in [Z]).
This value V is about 18.94... and is obtained for z = 4.52299.. (twentieth zero).
It appears that H(4.52299..) is negative.
Thus ff reaches a maximum if xx[1] = ... = xx[1000] = 4.52299..  and then takes the value 
C + 18.94.. * 1000 (which is the highest value of ff within [-5, 5]^1000

If I'm not mistaken

Here is an example (data XY come from the CurveFitting:-Spline help page)
The trap would be to search the extrema of each pice of the global function without accounting for the interval this piece is defined on. Another trap would be to search for extrema outsise the interval defined by the data.

 

restart:

with(CurveFitting):

XY := [[0,0],[1,1],[2,4],[3,3]];
S := Spline([[0,0],[1,1],[2,4],[3,3]], v);

XY := [[0, 0], [1, 1], [2, 4], [3, 3]]

 

piecewise(v < 1, (4/5)*v^3+(1/5)*v, v < 2, -2*v^3+(42/5)*v^2-(41/5)*v+14/5, (6/5)*v^3-(54/5)*v^2+(151/5)*v-114/5)

(1)

# ranges

bounds := min(op~(1, XY)), seq(rhs(op(k, S)), k in [seq(1..nops(S)-1, 2)]), max(op~(1, XY));
ranges := seq(bounds[k]..bounds[k+1], k=1..numelems([bounds])-1)

0, 1, 2, 3

 

0 .. 1, 1 .. 2, 2 .. 3

(2)

# 1st derivative of the S

dS := seq(diff(op(k, S), v), k in [seq(2..nops(S), 2)]), diff(op(-1, S), v);

(12/5)*v^2+1/5, -6*v^2+(84/5)*v-41/5, (18/5)*v^2-(108/5)*v+151/5

(3)

# extrema

sol := seq([fsolve(dS[k], v=ranges[k])], k=1..numelems([dS]))

[], [], [2.218264040]

(4)

# reconstruction as a piecewise solution

piecewise(seq(op([v >= op(1, ranges[k]) and v < op(2, ranges[k]), sol[k]]), k=1..numelems([sol])))

piecewise(0 <= v and v < 1, [], 1 <= v and v < 2, [], 2 <= v and v < 3, [2.218264040])

(5)

 


 

Download extrema.mw

@JAMET  @Kitonum

Hi JAMET, 

I don't know if the solution Kitonum has provided suits you (I understood you were looking for something more formal?).
If I'm mistaken forget what comes below, otherwise this worksheet might interest you.
From a line L defined by its cartesian equation it contains the step by step procedure to:

  1. construct the equation of L in the complex plane
  2. construct the equation of the line orthogonal to L and which passes through (0, 0)
  3. construct the expression of the intersection of these two lines

All of this is done formally, two examples follow with drawings

 

restart:

alias(conj = conjugate):

assumptions := (a::real, b::real, c::real, x::real, y::real)

a::real, b::real, c::real, x::real, y::real

(1)

# Cartesian equation of line L

L := a*x + b*y - c

a*x+b*y-c

(2)

# let Z a complex number

rel_1 := Z = x + I*y

Z = x+I*y

(3)

# let Z bar its conjugate

rel_2 := `#mrow(mover(mo(Z),mo("&#x305;")))` = eval(conj(Z), rel_1) assuming assumptions

`#mrow(mover(mo(Z),mo("&#x305;")))` = x-I*y

(4)

# explain x and y in terms of Z and Z bar

interface(warnlevel=0):
rel_3 := op( solve([rel_1, rel_2], [x, y]) )

[x = (1/2)*Z+(1/2)*`#mrow(mover(mo(Z),mo("&#x305;")))`, y = -((1/2)*I)*(Z-`#mrow(mover(mo(Z),mo("&#x305;")))`)]

(5)

# rewrite L by using the above equalities

eval(L, rel_3);

a*((1/2)*Z+(1/2)*`#mrow(mover(mo(Z),mo("&#x305;")))`)-((1/2)*I)*b*(Z-`#mrow(mover(mo(Z),mo("&#x305;")))`)-c

(6)

# take the numerator of this expression

numer(%);

-I*b*Z+I*b*`#mrow(mover(mo(Z),mo("&#x305;")))`+a*Z+a*`#mrow(mover(mo(Z),mo("&#x305;")))`-2*c

(7)

# collect the result according to Z and Z bar

`#mrow(mo("&#8466;"))`  := collect(%, [Z, `#mrow(mover(mo(Z),mo("&#x305;")))`])

(-I*b+a)*Z+(I*b+a)*`#mrow(mover(mo(Z),mo("&#x305;")))`-2*c

(8)

# define the complex number v from a and b

rel_4 := v = a+I*b

v = I*b+a

(9)

# let v bar its conjugate

rel_5 := `#mrow(mover(mo(v),mo("&#x305;")))` = eval(conj(v), rel_4)  assuming assumptions;

`#mrow(mover(mo(v),mo("&#x305;")))` = -I*b+a

(10)

# reverse rel_4 and rel_5

rel_6 := (rhs=~lhs)~([rel_4, rel_5])

[I*b+a = v, -I*b+a = `#mrow(mover(mo(v),mo("&#x305;")))`]

(11)

# here is the equation of L in the complex plane

`#mrow(mo("&#8466;"))` := eval(`#mrow(mo("&#8466;"))`, rel_6)

Z*`#mrow(mover(mo(v),mo("&#x305;")))`+v*`#mrow(mover(mo(Z),mo("&#x305;")))`-2*c

(12)

# find the complex equation of the line orthogonal to L and passing through the origin [0, 0]
#
# firstly define w as a rotation by Pi/2 of the vector whose affix is [a, b] in the complex plane

rel_7 := w = v*I

w = I*v

(13)

# let w bar its conjugate

rel_8 := `#mrow(mover(mo(w),mo("&#x305;")))` = eval(conj(w), rel_7);
rel_8 := eval(rel_8, conj(v) = `#mrow(mover(mo(v),mo("&#x305;")))`);

`#mrow(mover(mo(w),mo("&#x305;")))` = -I*conj(v)

 

`#mrow(mover(mo(w),mo("&#x305;")))` = -I*`#mrow(mover(mo(v),mo("&#x305;")))`

(14)

# here is the equation of the line orthogonal to L, and which passes through the origin

`#mrow(mo("&#8459;"))` := eval(
                                `#mrow(mo("&#8466;"))`,
                                [
                                  v=w,
                                  `#mrow(mover(mo(v),mo("&#x305;")))` = `#mrow(mover(mo(w),mo("&#x305;")))`,
                                  c=0
                                ]
                              )

Z*`#mrow(mover(mo(w),mo("&#x305;")))`+w*`#mrow(mover(mo(Z),mo("&#x305;")))`

(15)

# where do lines L and H intersect ?
#
# let's form this system in Z and Z bar (brute force)

sys := [ `#mrow(mo("&#8466;"))` = 0, eval(`#mrow(mo("&#8459;"))`, [rel_7, rel_8]) = 0 ]:

print~(sys):

Z*`#mrow(mover(mo(v),mo("&#x305;")))`+v*`#mrow(mover(mo(Z),mo("&#x305;")))`-2*c = 0

 

-I*Z*`#mrow(mover(mo(v),mo("&#x305;")))`+I*v*`#mrow(mover(mo(Z),mo("&#x305;")))` = 0

(16)

# where do lines L and H intersect ?

sol := op( solve( sys, [ Z, `#mrow(mover(mo(Z),mo("&#x305;")))` ] ) )

[Z = c/`#mrow(mover(mo(v),mo("&#x305;")))`, `#mrow(mover(mo(Z),mo("&#x305;")))` = c/v]

(17)

# use the expression of v to explicit the solution (obviously Z and Z bar are conjugate)

sol := eval(sol, [rel_4, rel_5]);

[Z = c/(-I*b+a), `#mrow(mover(mo(Z),mo("&#x305;")))` = c/(I*b+a)]

(18)

# transform the solution into a prettier form

sol := lhs~(sol) =~ numer~(rhs~(sol)) *~ conj~(denom~(rhs~(sol)))
                    /~
                    expand~(denom~(rhs~(sol)) *~ conj~(denom~(rhs~(sol))) )  assuming assumptions

[Z = -c*(-I*b-a)/(a^2+b^2), `#mrow(mover(mo(Z),mo("&#x305;")))` = c*(-I*b+a)/(a^2+b^2)]

(19)

# this complex number -represents the affix of the intersection point of the two lines

`#mrow(mo("&#8484;"))` := simplify(eval(Z, sol));  # asuming a^2+b^2 <> 0

c*(I*b+a)/(a^2+b^2)

(20)

# example

data := [a=3, b=-2, c=1, x=X, y=Y]:

with(plots):

f := eval(a*X + b*Y - c, data);
`#mrow(msub(mo(f),mo("&#8869;")))` := simplify(
                                        eval(
                                          eval(
                                            eval(`#mrow(mo("&#8459;"))`, [rel_1, rel_2, rel_7, rel_8]),
                                            [rel_4, rel_5]
                                          ),
                                          data
                                        )
                                      );
__Z := eval(`#mrow(mo("&#8484;"))`, data);

display(
  implicitplot(f, X=-1..1, Y=-1..1, color=blue),
  implicitplot(`#mrow(msub(mo(f),mo("&#8869;")))`, X=-1..1, Y=-1..1, color=red),
  pointplot([[Re(__Z), Im(__Z)]], symbol=circle, symbolsize=20),
  gridlines=true,
  scaling=constrained
);

3*X-2*Y-1

 

6*Y+4*X

 

3/13-(2/13)*I

 

 

# just another example

data := [a=0, b=-3, c=-2, x=X, y=Y]:

with(plots):

f := eval(a*X + b*Y - c, data);
`#mrow(msub(mo(f),mo("&#8869;")))` := simplify(
                                        eval(
                                          eval(
                                            eval(`#mrow(mo("&#8459;"))`, [rel_1, rel_2, rel_7, rel_8]),
                                            [rel_4, rel_5]
                                          ),
                                          data
                                        )
                                      );
__Z := eval(`#mrow(mo("&#8484;"))`, data);

display(
  implicitplot(f, X=-1..1, Y=-1..1, color=blue),
  implicitplot(`#mrow(msub(mo(f),mo("&#8869;")))`, X=-1..1, Y=-1..1, color=red),
  pointplot([[Re(__Z), Im(__Z)]], symbol=circle, symbolsize=20),
  gridlines=true,
  scaling=constrained
);

-3*Y+2

 

6*X

 

(2/3)*I

 

 

 


 

Download formal_solution_in_complex_plane.mw

 

Hi, 
Have you try to use the geometry package?
Here are some of its capabilities that could interest you.


 

restart:

with(geometry):
with(RealDomain):

local D, O:

_EnvHorizontalName := x: _EnvVerticalName := y:
point(O, [0, 0]);

O

(1)

# D is assumed neither horizontal nor vertical

line(D, a*x+b*y+c) assuming a^2 + b^2 > 0, a<>0, b<>0:
PerpendicularLine(Perp, O, D):
projection(H, O, D):
coordinates(H);

[-a*c/(a^2+b^2), -b*c/(a^2+b^2)]

(2)

# D is assumed vertical

line(D, b*y+c) assuming b<>0, b^2 > 0:
PerpendicularLine(Perp, O, D):
projection(H, O, D):
coordinates(H);

[0, -c/b]

(3)

# D is assumed horizontal

line(D, a*x+c) assuming a<>0, a^2 > 0:
PerpendicularLine(Perp, O, D):
projection(H, O, D):
coordinates(H);

[-c/a, 0]

(4)

 


 

Download Pied_de_nez.mw

 

There are probably a lot of ways to do that.
Here is one of them
(the set of points I construct are in some sense the "shortest" for adding them any other point from the initial list obviously creates a new set of points whose sum id larger than the target sum)

 

restart:

randomize():

M    := 10.:
roll := rand(0. .. M):
N    := 10:
L    := [seq(roll(), k=1..N)]

[3.829591035, 3.239717830, 8.629534219, 7.550798661, 7.992442854, 7.267252468, 2.715869968, 3.288246716, 3.345034572, 4.004192509]

(1)

ROLL := rand(0..N*M):
S    := ROLL();  # target sum

8.585001451

(2)

# First way

c := Statistics:-CumulativeSum(L):
if c[-1] < S then
  printf("The sum of all the numbers (%a) is less than rge desired sum (%a)\n", c[-1], S)
else
  select[flatten](`<`, c, S);
  selected_numbers := L[1..numelems(%)+1];

  # verify
  is(add(selected_numbers) > S)
end if;

Array([3.82959103500000, 7.06930886500000])

 

[3.829591035, 3.239717830, 8.629534219]

 

true

(3)

# Repeat this many (R) times.

R := 10:
for r from 1 to R do
  rL := combinat:-randperm(L):
  c := Statistics:-CumulativeSum(rL):
  if c[-1] < S then
    printf("The sum of all the numbers (%a) is less than rge desired sum (%a)\n", c[-1], S)
  else
    select[flatten](`<`, c, S);
    selected_numbers := rL[1..numelems(%)+1];
    printf("%a\n", selected_numbers);
  end if;
end do:

[2.715869968, 8.629534219]
[4.004192509, 7.267252468]
[8.629534219]
[3.829591035, 3.288246716, 8.629534219]
[2.715869968, 3.829591035, 7.267252468]
[8.629534219]
[7.550798661, 7.992442854]
[3.829591035, 7.267252468]
[7.267252468, 7.992442854]
[8.629534219]

 

 

 

 

 

 

# select the smallest number in wich is larger than 3.14

L[ListTools:-BinaryPlace(L, 3.14)+1];
min(select[flatten](`>`, L, 3.14))

4

 

4

(4)

# select the smallest number in wich is larger or equal than 4

min(select[flatten](`>`, L, 4));
min(select[flatten](`>=`, L, 4));

5

 

4

(5)

# all the points in L that are larger or equal than 4


LL := select[flatten](`>=`, L, 4)
 

[4, 5, 6, 7, 8, 9, 10]

(6)

# randomly select some of them, let's say 3 of them

combinat:-randperm(LL)[1..3]

[5, 4, 10]

(7)

 

 

 


 

Download search.mw

Without any constraints on the selcted numbers, all the solutions are equivalent.
But maybe you are looking for the selected numbers whose sum is the closest to the target sum?
Or the minimum set of numbers whose sum is larger than this target sum?

 

interface(imaginaryunit=i)

A minimal answer (LinearFit can return a lot of informations, please look to the corresponding help page).

Download LinearFit.mw

Probably the shortest coding (note that prime number 2 is congruent to 2, so there also exists a list L[2] with only one member)

Statistics:-Tally(`mod`(select(isprime, L), 8))

If you do not want to compute the number of elements L[2] has (which is obviously 1 whatever the length of L), the simplest thing do do is to begin the construction from the 3rd element of L (the first (1) is not prime, and the second (2), is the prime you don't want to count)

Statistics:-Tally(`mod`(select(isprime, L[3..-1]), 8))
 

 

@minhhieuh2003 

I have no cure to that, only an explanation.
If you type whattype(1/x) you will get `^` meaning 1/x is not interpreted as a fraction but as somthing (x) raised to some power (-1).
This is confirmed by type(1/x, fraction) whose output is  false.
(but type(2/3, fraction)returns  true).
Then the latex command probably uses the representation of 1/x to do the LaTeX translation, and this representation is of power type and not of fraction type.
More precisely I think latex uses the representation produced by the command dismantle :
dismantle(1/x)
PROD(3)
   NAME(4): x
   INTNEG(2): -1

You see here that 1/x is some name raised to an integer power equal to -1. This is consistent with the outputs of whattype(1/x) and type(1/x, fraction).

Look now what dismantle(2/3) returns
dismantle(2/3)
RATIONAL(3): 2/3
   INTPOS(2): 2
   INTPOS(2): 3


To force latex to return something like \frac{1}{x} you could think to write your own procedure, for instance 
`latex/frac` := proc()
  cat(`\\frac{`, `latex/print`(args[1]), `}{`, `latex/print`(args[2]), `}` )
end proc;

See what the result will be
printf("%s", `latex/frac`(1, x^2+2))

But the harder problem will be to force maple to interpret  1/expression as a fraction and not as this expression to the power -1 (more generally to interpret expression_1/expression_2 as a fraction).
And this one I'm not sure it can be solved ???


As a last resort, you can always write your own Maple --> LateX converter.
You fill find in the attached file ad hoc a snippet of an ad hoc LaTeX translator 
 

restart

# Ultimately you can do that by hand:

cat(`\\int{2}{3}\\frac{1}{`, `latex/print`(x^2+2), `}=`, `latex/print`(eval(parse(expr))) ):
printf("%s", %);

\int{2}{3}\frac{1}{{x}^{2}+2}={\it expr}

 

# But it would be more generic to write your own procedure
 

`latex/frac` := proc()
  if patmatch(args[1], (a::anything)/(b::anything), 's' ) then
    if rhs(s[2]) <> 1 then
      cat(`\\frac{`, `latex/print`(rhs(s[1])), `}{`, `latex/print`(rhs(s[2])), `}` )
    else
      cat(`latex/print`(rhs(s[1])))
    end if:
  end if:
end proc:

printf("%s", `latex/frac`(1/(x^2+2)))

\frac{1}{{x}^{2}+2}

 

printf("%s", `latex/frac`(x*(x^2+2)))

x \left( {x}^{2}+2 \right)

 

 

 

 

 

`LateX/Int` := proc(expr)
  local L, i, v, lob, upb:
  L := [op(0..-1, expr)]:
  if L[1] <> Int then error "the expression is not of type Int" end if:
  if patmatch(L[2], (a::anything)/(b::anything), 's' ) then
    if rhs(s[2]) <> 1 then
      i := cat(`\\frac{`, `latex/print`(rhs(s[1])), `}{`, `latex/print`(rhs(s[2])), `}` )
    else
      i := cat(`latex/print`(rhs(s[1])))
    end if:
  end if:
  v := lhs(L[3]):
  if type(L[3], `=`) then
    lob := op(1, rhs(L[3])):
    upb := op(2, rhs(L[3])):
    cat(`\\int_{`, `latex/print`(lob), `}^{`, `latex/print`(upb), `}\\!`, i, `\,{\\rm d}`, `latex/print`(v))
  else
    cat(`\\int_{`, `latex/print`(eval(lob)), `}^{`, `latex/print`(upb), `}\\!`, i, `\,{\\rm d}`, `latex/print`(v))
  end if:
end proc:

latex(Int(x, x=0..1))

\int_{0}^{1}\!x\,{\rm d}x

 

printf("%s", `LateX/Int`(Int(x, x=0..1)))

\int_{0}^{1}\!x,{\rm d}x

 

printf("%s", `LateX/Int`(Int(1/(x^2+2), x=0..1)))

\int_{0}^{1}\!\frac{1}{{x}^{2}+2},{\rm d}x

 

 


 

Download latex.mw

 

 

Hi, 


I guess that R1 and R2 have bounded supports S1 and S2 (if not either the min or the max or both of them is infinite).
I also admit that the function f : (R1, R2) ---> X is everywhere defined over S1xS2 (wich is not the case if f(R1,R2)=R1/R2 and S2 contains {0}).

So, if I understand correctly your question,  there is no need to use random variables to find the min and the max of X.

Look to the package Tolerances, it should answer your question.

The case of "nom" (I guess "nominal") is more subtle. But if you do not rely "nominal" to some statististics of R1 and R2 ("nominal" is just some specific value, not a mean, median, mode, ...) then Tolerances is stiil the paxkage to use.

If "nominal" repesent mean values (for instance), the "nominal" value of X can be easily obtained by using the package ScientificErrorAnalysis.
Feel free to tell me if this didn't help

If you want to represent a univariate polynomial in some orthogonal polynomial expansion, this might interest you.

PS: the last command InsertContent(Worksheet(Group(Input( A ))));  can't be displayed here, just run the code to visualize it.

restart:

interface(version);

`Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895`

(1)

with(DocumentTools):
with(DocumentTools:-Layout):

OP := with(orthopoly)

[G, H, L, P, T, U]

(2)

F := proc(P, x, B)
  local d, C, Q, symb:

  if _npassed = 4 and _passed[3] <> G then
    error cat("A fourth argument (", _passed[4], ") is provided but the orth. poly is not of type G")
  end if;
  if _npassed = 3 and _passed[3] = G then
    error cat("Using type G orth. poly. needs to pass F a fourth argument, see orthopoly[G] help page")
  end if;
  d    := degree(P, x):
  symb := convert(cat("#mrow(mo(", B, "))"), name);
  if _passed[3] <> G then
    Q := add(a[n]*B(n, x), n=0..d):
    C := op(solve( [ coeffs( collect(expand(P-Q), x), x) ], [ seq(a[n], n=0..d) ] )):
    return add(rhs~(C) *~ [ seq(symb(n, x), n=0..d) ])
  else
    Q := add(a[n]*B(n, _passed[4], x), n=0..d):
    C := op(solve( [ coeffs( collect(expand(P-Q), x), x) ], [ seq(a[n], n=0..d) ] )):
    return add(rhs~(C) *~ [ seq(symb(n, _passed[4], x), n=0..d) ])
  end if:
end proc:

# Laguerre representation of polynom Pol
# (watchout : do not name P this polunomial to avoid conflicts with orthopoly P)
# note the character "L" that appears here is not the "L" which represents LaguerreL polynomials

Pol := randpoly(x, dense, degree = 4);

B   := L:
pol := F(Pol, x, B);
pol := F(Pol, x, B, 1);

print();
B   := G:
pol := F(Pol, x, B);
pol := F(Pol, x, B, 1);

-83*x^4+98*x^3-48*x^2-19*x+62

 

-1457*`#mrow(mo(L))`(0, x)+6415*`#mrow(mo(L))`(1, x)-10284*`#mrow(mo(L))`(2, x)+7380*`#mrow(mo(L))`(3, x)-1992*`#mrow(mo(L))`(4, x)

 

Error, (in F) A fourth argument (1) is provided but the orth. poly is not of type G

 

 

Error, (in F) Using type G orth. poly. needs to pass F a fourth argument, see orthopoly[G] help page

 

(317/8)*`#mrow(mo(G))`(0, 1, x)+15*`#mrow(mo(G))`(1, 1, x)-(441/16)*`#mrow(mo(G))`(2, 1, x)+(49/4)*`#mrow(mo(G))`(3, 1, x)-(83/16)*`#mrow(mo(G))`(4, 1, x)

(3)

 

# Extract the symbol used and replace it by B and verify the result is Pol

symb := op(0, indets(pol, function)[1]);
eval(pol, symb=B);
%-Pol

`#mrow(mo(G))`

 

-83*x^4+98*x^3-48*x^2-19*x+62

 

0

(4)

# Several orthogonal polynomials representations
# Watchout: Gegenbauer polynomials ar not handled
N := numelems(OP):
M := Matrix(N+1, 2):
M[1, 1] := 'canonical':
M[1, 2] := Pol:
k := 2:
for B in OP do
  M[k, 1] := B:
  M[k, 2] := `if`(B <> G, F(Pol, x, B), F(Pol, x, B, 1)); # for instance
  k := k+1:
end do:


A := Table(
            seq(Column(), k=1..12), widthmode=percentage, width=60, # width has to be adjusted for larger expressions
            seq(
                 Row(
                      Cell( Textfield( style=TwoDimOutput, Equation(M[k,1]) ), columnspan=2),
                      Cell( Textfield( style=TwoDimOutput, Equation(M[k, 2]) ), columnspan=10)
                 ),
                 k=1..N+1
            )
     ):
InsertContent(Worksheet(Group(Input( A ))));

  

 


 

Download OrthopolyRepresentation.mw


You will find below how to modify the structure plotted by DrawGraph in order that the thickness of the edges are proportional to their weights.
The code is not robust because it doesn't gandle situations where the weights are not integers nor situations where the range of the weight is very large (let's say 1..100 for instance).
I tried to be pedagocic by diplaying intermediate steps.

PS: I used MAPLE 2015... maybe more recent version have embedded features to do all thus stuff



 

restart:

interface(version);

`Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895`

(1)

with(GraphTheory):

G := Graph({[{1,2},2],[{2,3},1],[{1,3},3]});

GRAPHLN(undirected, weighted, [1, 2, 3], Array(%id = 18446744078219989110), `GRAPHLN/table/1`, Matrix(%id = 18446744078220584678))

(2)

p := DrawGraph(G):
plots:-display(p);

 

# What does p contain?

d := [op(p)]:
print~(d):

POLYGONS([[.5252013392, 1.025201339], [.5252013392, .9747986607], [.4747986606, .9747986607], [.4747986606, 1.025201339]], [[.9582140413, .2752013392], [.9582140413, .2247986606], [.9078113627, .2247986606], [.9078113627, .2752013392]], [[0.9218863736e-1, .2752013393], [0.9218863736e-1, .2247986607], [0.4178595884e-1, .2247986607], [0.4178595884e-1, .2752013393]], COLOR(RGB, 1, 1, .2, 1, 1, .2, 1, 1, .2), STYLE(PATCHNOGRID))

 

TEXT([.4999999999, 1.0], 1, FONT(HELVETICA, BOLD, 12))

 

TEXT([.9330127020, .2499999999], 2, FONT(HELVETICA, BOLD, 12))

 

TEXT([0.6698729810e-1, .2500000000], 3, FONT(HELVETICA, BOLD, 12))

 

POLYGONS([[.9330127020, .2499999999], [.4999999999, 1.0]], COLOR(RGB, 0, 0, 1), THICKNESS(2), STYLE(LINE))

 

POLYGONS([[0.6698729810e-1, .2500000000], [.4999999999, 1.0]], COLOR(RGB, 0, 0, 1), THICKNESS(2), STYLE(LINE))

 

POLYGONS([[0.6698729810e-1, .2500000000], [.9330127020, .2499999999]], COLOR(RGB, 0, 0, 1), THICKNESS(2), STYLE(LINE))

 

TEXT([.7613076212, .5508660253], "2", ALIGNRIGHT, ALIGNABOVE, FONT(HELVETICA, 11))

 

TEXT([.2469423789, .5461028857], "3", ALIGNRIGHT, ALIGNBELOW, FONT(HELVETICA, 11))

 

TEXT([.4133974597, .2543301270], "1", ALIGNRIGHT, ALIGNABOVE, FONT(HELVETICA, 11))

 

SCALING(CONSTRAINED)

 

AXESSTYLE(NONE)

(3)

nv := numelems(Vertices(G));
ne := numelems(Edges(G));

3

 

3

(4)

# d[1]                   draws the nv yellow squares that represent the vertices
# d[2]...d[nv+1]         draw the labels of the vertices
# d[nv+2]..d[nv+2+ne-1]  draw the ne edges
# ...
#
# q contains the elements of d that we won't change
# r contains the elements of d that we are going to change

q := d[[$1..nv+1, $nv+2+ne..numelems(d)]]:
print~(q):

print():
r := d[[$nv+2..nv+2+ne-1]]:
print~(r):

POLYGONS([[.5252013392, 1.025201339], [.5252013392, .9747986607], [.4747986606, .9747986607], [.4747986606, 1.025201339]], [[.9582140413, .2752013392], [.9582140413, .2247986606], [.9078113627, .2247986606], [.9078113627, .2752013392]], [[0.9218863736e-1, .2752013393], [0.9218863736e-1, .2247986607], [0.4178595884e-1, .2247986607], [0.4178595884e-1, .2752013393]], COLOR(RGB, 1, 1, .2, 1, 1, .2, 1, 1, .2), STYLE(PATCHNOGRID))

 

TEXT([.4999999999, 1.0], 1, FONT(HELVETICA, BOLD, 12))

 

TEXT([.9330127020, .2499999999], 2, FONT(HELVETICA, BOLD, 12))

 

TEXT([0.6698729810e-1, .2500000000], 3, FONT(HELVETICA, BOLD, 12))

 

TEXT([.7613076212, .5508660253], "2", ALIGNRIGHT, ALIGNABOVE, FONT(HELVETICA, 11))

 

TEXT([.2469423789, .5461028857], "3", ALIGNRIGHT, ALIGNBELOW, FONT(HELVETICA, 11))

 

TEXT([.4133974597, .2543301270], "1", ALIGNRIGHT, ALIGNABOVE, FONT(HELVETICA, 11))

 

SCALING(CONSTRAINED)

 

AXESSTYLE(NONE)

 

 

POLYGONS([[.9330127020, .2499999999], [.4999999999, 1.0]], COLOR(RGB, 0, 0, 1), THICKNESS(2), STYLE(LINE))

 

POLYGONS([[0.6698729810e-1, .2500000000], [.4999999999, 1.0]], COLOR(RGB, 0, 0, 1), THICKNESS(2), STYLE(LINE))

 

POLYGONS([[0.6698729810e-1, .2500000000], [.9330127020, .2499999999]], COLOR(RGB, 0, 0, 1), THICKNESS(2), STYLE(LINE))

(5)

# verification
plots:-display(r);

 

# get tertex positions


v := GetVertexPositions(G);

[[.4999999999, 1.0], [.9330127020, .2499999999], [0.6698729810e-1, .2500000000]]

(6)

# change the default thickness by the weight of the edge(assuming it is a positive number)

W := WeightMatrix(G);
S := NULL:

for s in r do
  __edge      := op(1, s):                           # the edge "s"
  __from_to   := map(ListTools:-Search, __edge, v):  # starts at "from" and ends at "to"
  __thickness := W[op(__from_to)]:
  
  S := S, subs(THICKNESS(2) = THICKNESS(__thickness), s)
end do:

PLOT(S);

W := Matrix(3, 3, { sparse_data }, storage = sparse, shape = [symmetric])

 

 

# plot the sequence q, S

PLOT(q[], S)

 

 


 

Download WeightedEdges.mw

 

To complete Kitonum's and Rouben's replies, and if you have Maple 2018 or higher, you can also use the Kriging method (from memory look to the CurveFitting package ... If I'm mistaken please corrrect me).
Less immediate if you do not know the basis of Kriging, but this will give you a more flexible way to obtain a function which interpolates your unknown function from a set of X-Y couples.
For I only have Maple 2015 right now, I've written an extremely simple procedure to do "simple Kriging". Here is a few examples of interpolators you can build (the X-Y couples come from your plot)?
They differ by the value of a scalar parameter which controls the smoothness of the interpolation.

@tayyild 

Here is a solution which works for 2 or more polynomials of arbitrary number of monomials.
The procedure needs two lists P and V of equal lengths:

  • P[n] denotes the number of monomials of polynom n
  • V[n] denotes the symbol used to represent polynom n

Your first example is coded P := [2, 2], V=[x, y] and your second P := [3, 3], V=[x, y].

A_general_solution.mw

First 48 49 50 51 52 53 54 Last Page 50 of 62