Advanced Search

vv

13992 Reputation

20 Badges

10 years, 38 days
Contact vv

MaplePrimes Activity


These are answers submitted by vv

I don't think it has a special name;...

vv 13992
January 26 2018
0 2

I don't think it has a special name; it is a multiplication by a diagonal matrix.

But are you sure it is useful? AFAIK the eigenvector algorithms always do a row (or column) normalization.

evalc...

vv 13992
January 25 2018
2 1

Use
evalc(Re(expr));

 

Solution...

vv 13992
January 25 2018
3 3
> 

restart;
> 

eq1 := diff(u1(x), x, x)+diff(u2(x), x)+int(2*x*s*(u1(s)-3*u2(s)), s = 0 .. 1) = 6*x^2+3*x*(1/10)+8;
eq2 := diff(u1(x), x)+diff(u2(x), x, x)+int((3*(s^2+2*x))*(u1(s)-2*u2(s)), s = 0 .. 1) = 21*x+4/5;
bcs := u1(0)+(D(u1))(0) = 1, u2(0)+(D(u2))(0) = 1, u1(1)+(D(u1))(1) = 10, u2(1)+(D(u2))(1) = 7;

EQ1 := diff(u1(x), x, x)+diff(u2(x), x)+a*x = 6*x^2+3*x*(1/10)+8;
EQ2 := diff(u1(x), x)+diff(u2(x), x, x)+b*x+c = 21*x+4/5;

diff(diff(u1(x), x), x)+diff(u2(x), x)+int(2*x*s*(u1(s)-3*u2(s)), s = 0 .. 1) = 6*x^2+(3/10)*x+8

  

diff(u1(x), x)+diff(diff(u2(x), x), x)+int(3*(s^2+2*x)*(u1(s)-2*u2(s)), s = 0 .. 1) = 21*x+4/5

  

u1(0)+(D(u1))(0) = 1, u2(0)+(D(u2))(0) = 1, u1(1)+(D(u1))(1) = 10, u2(1)+(D(u2))(1) = 7

  

diff(diff(u1(x), x), x)+diff(u2(x), x)+a*x = 6*x^2+(3/10)*x+8

  

diff(u1(x), x)+diff(diff(u2(x), x), x)+b*x+c = 21*x+4/5

 (1)
> 

sol:=dsolve({EQ1,EQ2},{u1(x),u2(x)}):
> 

U1:=unapply( eval(u1(x),sol),x):
U2:=unapply( eval(u2(x),sol),x):
> 

sys:=eval([(lhs-rhs)(eq1),(lhs-rhs)(eq2),bcs],[u1=U1,u2=U2]):
> 

map(coeffs,sys,x):
> 

consts:=solve(%):
> 

SOL:=simplify(eval(sol,consts));

{u1(x) = (-55897*exp(1-x)+27844*exp(2-x)+35538*exp(1+x)+(229194*x^2+1583820*x-1422696)*exp(1)+(-45192*x^2-385320*x+335340)*exp(2)-300822*x^2-1458780*x-55263*exp(x)+1319226)/(-49980*exp(2)+232200*exp(1)-250080), u2(x) = (-55897*exp(1-x)+27844*exp(2-x)-35538*exp(1+x)+(464400*x^3-681738*x^2+1399212*x-1095936)*exp(1)+(-99960*x^3+160164*x^2-309456*x+259476)*exp(2)-500160*x^3+669894*x^2-1398996*x+55263*exp(x)+1038390)/(-49980*exp(2)+232200*exp(1)-250080)}

 (2)
> 

U1:=unapply( eval(u1(x),SOL),x):             # Check
U2:=unapply( eval(u2(x),SOL),x):             #
simplify(eval({eq1,eq2,bcs},[u1=U1,u2=U2])); #

{1 = 1, 7 = 7, 10 = 10, 21*x+4/5 = 21*x+4/5, 6*x^2+(3/10)*x+8 = 6*x^2+(3/10)*x+8}

 (3)
> 

 


Download sysode-int.mw

By recurrence...

vv 13992
January 24 2018
1 1

 

n is supposed to be a positive integer.

We compute the indefinite integral.

 

 

> 

J:=n -> Int((2*cos(Pi*x*n/T)^2-1)*sin(Pi*x*n/T)*cos(Pi*x*n/T)/(sin(Pi*x/T)*cos(Pi*x/T)), x);

proc (n) options operator, arrow; Int((2*cos(Pi*x*n/T)^2-1)*sin(Pi*x*n/T)*cos(Pi*x*n/T)/(sin(Pi*x/T)*cos(Pi*x/T)), x) end proc

 (1)
> 

value(combine(J(n+1)-J(n)));

(1/2)*T*sin(2*Pi*x*(2*n+1)/T)/(Pi*(2*n+1))

 (2)
> 

K:=unapply(%,n);

proc (n) options operator, arrow; (1/2)*T*sin(2*Pi*x*(2*n+1)/T)/(Pi*(2*n+1)) end proc

 (3)
> 

J(n)=value(J(1))+Sum(K(k),k=1..n-1);  # Answer

Int((2*cos(Pi*x*n/T)^2-1)*sin(Pi*x*n/T)*cos(Pi*x*n/T)/(sin(Pi*x/T)*cos(Pi*x/T)), x) = -x+2*T*((1/2)*cos(Pi*x/T)*sin(Pi*x/T)+(1/2)*Pi*x/T)/Pi+Sum((1/2)*T*sin(2*Pi*x*(2*k+1)/T)/(Pi*(2*k+1)), k = 1 .. n-1)

 (4)
> 

value(%);  # Version of the answer

int((2*cos(Pi*x*n/T)^2-1)*sin(Pi*x*n/T)*cos(Pi*x*n/T)/(sin(Pi*x/T)*cos(Pi*x/T)), x) = -x+2*T*((1/2)*cos(Pi*x/T)*sin(Pi*x/T)+(1/2)*Pi*x/T)/Pi+((1/8)*I)*T*(exp((2*I)*Pi*x*(2*n+1)/T)*LerchPhi(exp((4*I)*Pi*x/T), 1, n+1/2)-exp((6*I)*Pi*x/T)*LerchPhi(exp((4*I)*Pi*x/T), 1, 3/2)+exp(-(6*I)*Pi*x/T)*LerchPhi(exp(-(4*I)*Pi*x/T), 1, 3/2)-LerchPhi(exp(-(4*I)*Pi*x/T), 1, n+1/2)*exp(-(2*I)*Pi*x*(2*n+1)/T))/Pi

 (5)
> 

simplify(diff( (lhs-rhs)(%), x)); # Check

0

 (6)

Maths...

vv 13992
January 23 2018
0 0

This is a maths problem, not a Maple one. See https://en.wikipedia.org/wiki/Cycles_and_fixed_points

Of course it is easy to implement any of the formulae.

Direct solution...

vv 13992
January 23 2018
2 7

Here is a direct solution (based on simplex) asked by the OP who has an older version of Maple.

FacetsOfPolyhedralCone:=proc(L::list)
local i,j,J,A,X,u;
X:=indets(L,name);
A:=LinearAlgebra:-GenerateMatrix(L,X)[1]:
J:={seq(1..nops(L))}:
for j in J do
  add(A[j]-u[i]*A[i],i=J minus {j}); convert(%,list)=~0;
  if simplex[minimize](0, %, NONNEGATIVE) <>{} then J:=J minus {j} fi;
od:
J, L[[J[]]];
end:

L := [y-z, 3*y-2*z, 2*y-2*z, x-2*y+z, x-y, 2*x-y, x-z, x+y-z, x, y, z]:  # >= 0

FacetsOfPolyhedralCone(L);
       {3, 4, 11}, [2*y-2*z, x-2*y+z, z]

Bug+...

vv 13992
January 22 2018
0 1 
solve({x[5] = x[2]/x[1], x[6] = x[3]/x[2], x[7] = x[1]/x[4],
       x[8] = (2*x[2]+x[4])/(2*x[1]+x[3]+x[4])}, {x[1], x[2], x[3], x[4]});

         {x[1] = 0, x[2] = 0, x[3] = 0, x[4] = 0}

It is a bug. The result of eliminate sho...

vv 13992
January 22 2018
1 0

It is a bug. The result of eliminate should be:

Mathematically it should be also added: x[5]<>0, x[7]<>0, x[4]<>0,  x[6] <> -(2*x[7]+1)/(x[5]*x[7])

@mmcdara  Let's do the computa...

vv 13992
January 21 2018
0 4

@mmcdara 

Let's do the computations for 3 digits.

restart;
Digits:=3:
a,b,c:= sqrt~([2,3,6])[]:
fa,fb,fc := evalf(%);
                 fa, fb, fc := 1.41, 1.73, 2.45
fa*fb;
                              2.44
# Do you expect  2.45?  How?
141*173;
                             24393

 

 

?...

vv 13992
January 21 2018
0 0 
num1dsubspaces := (p,n) -> (p^n - 1)/(p-1)

 

PolyhedralSets...

vv 13992
January 20 2018
2 3 
ineqs:=
y-z >=0,
3*y-2*z >=0,
2*y-2*z >=0,
x-2*y+z >=0,
x-y >=0,
2*x-y >=0,
x-z >=0,
x+y-z >=0,
x>=0,
y>=0,
z>=0:
with(PolyhedralSets):
A := PolyhedralSet([ineqs]);
B := PolyhedralSet([ineqs[3..4],ineqs[-3..-1]]);


evalb(convert(A,string)=convert(B,string));

        true

plot (smaller interval)...

vv 13992
January 20 2018
0 3

Add:

f1,f2,f3:=seq(unapply('evalf'(x__||i),t2), i=1..3):
plot([f1,f2,f3], 0..0.2);

You will not be able to plot over  0 .. Pi/2  because the integrals are not convergent here.

I think that for a listlist L, only the ...

vv 13992
January 20 2018
0 3

I think that for a listlist L, only the following constructs L[u,v] should be accepted

L[i1..i2, j1..j2],  L[i, j1..j2],  L[i1..i2, j],  L[i, j]

(they are consistent with rtables).

The other ones such as L[[i1,i2,...], [j1,j2,...]] should produce exceptions (errors)
because they are not consistent with rtables and they are not needed: L[u,v]  reduces to L[u][v] in these cases.
 

compiled - fixed...

vv 13992
January 20 2018
1 1


 

> 

MinVal := proc(p :: Array(datatype=integer[4])
               , lg :: Array(datatype=float[8])
               , pmin :: Array(datatype=integer[4])
               , m :: float
               , K :: float
               , n :: posint
               , len :: posint
              )
local i,j,v,x;
    v := 0.0;
    for i from 0 to n-1 by len do
        x := 0.0;
        for j to len do
            x := x + lg[p[i+j]];
        end do;
        v := v + abs(x - K);
    end do;

    # test for minimum
    if m <= v then
        v := m;
    else
        # update pmin
        for i to n do
            pmin[i] := p[i];
        end do;
    end if;
    return v;
end proc:
> 

Var := proc(L::list(positive)
            , len :: posint
            , num :: posint
           )
local M, P, i, j, K, lg, n, p, pmin;
    n := numelems(L);
    if len*num <> n then
        error "invalid partition specification";
    end if;

    # allocate arrays used by MinVal
    pmin := Array(1..n, 'datatype'=integer[4]);

    lg := Array(evalf(ln~(L)),'datatype'=float[8]);
    K := add(i, i=lg)/num;
    M := K*len;

    P := Iterator:-SetPartitionFixedSize([len$num]);

    for p in P do
        M := MinVal(p, lg, pmin, M, K, n, len);
    end do;

    return ListTools:-LengthSplit(convert(pmin,list), len);

end proc:
> 

L := [1829.0, 1644.0, 1594.0, 1576.0, 1520.0, 1477.0, 1477.00, 1404.0
      , 1392.0, 1325.0, 1313.0, 1297.0, 1292.0, 1277.0, 1249.0, 1236.0]:
 
> 

MinVal := Compiler:-Compile(MinVal):
CodeTools:-Usage((Var)(L, 4, 4));

memory used=405.19MiB, alloc change=-1.00MiB, cpu time=8.58s, real time=8.59s, gc time=109.20ms

  

[1, 9, 13, 15], [2, 4, 14, 16], [3, 6, 10, 11], [5, 7, 8, 12]

 (1)
> 

indxs:=[%[1..4]]:
vals := [A||(1..16)]:
map(i->vals[i], indxs);
map(i->L[i], indxs);

[[A1, A9, A13, A15], [A2, A4, A14, A16], [A3, A6, A10, A11], [A5, A7, A8, A12]]

  

[[1829.0, 1392.0, 1292.0, 1249.0], [1644.0, 1576.0, 1277.0, 1236.0], [1594.0, 1477.0, 1325.0, 1313.0], [1520.0, 1477.00, 1404.0, 1297.0]]

 (2)
> 

 


 

Download compiled.mw

plot3d...

vv 13992
January 19 2018
1 1

plot3d wants a list of 3 procedures but receives a single procedure (which is list-valued).

First 77 78 79 80 81 82 83 Last Page 79 of 120