Advanced Search

vv

14107 Reputation

20 Badges

10 years, 132 days
Contact vv

MaplePrimes Activity


These are answers submitted by vv

scalars...

vv 14107
June 05 2017
2 3


 

It seems that now the Vector and Matrix constructors accept only scalars.

But we may use delayed evaluation.

 

> 

Vector(2,(a) -> 'Matrix'(2,2,(b,c) -> m||a||b||c));

Vector(2, {(1) = Matrix(2, 2, {(1, 1) = m111, (1, 2) = m112, (2, 1) = m121, (2, 2) = m122}), (2) = Matrix(2, 2, {(1, 1) = m211, (1, 2) = m212, (2, 1) = m221, (2, 2) = m222})})

 (1)
> 

A:=<1,2;3,4>;B:=<2,3,4>;

A := Matrix(2, 2, {(1, 1) = 1, (1, 2) = 2, (2, 1) = 3, (2, 2) = 4})

  

Vector[column](%id = 18446744074182844158)

 (2)
> 

Vector([B,B]); #assemble

Vector(6, {(1) = 2, (2) = 3, (3) = 4, (4) = 2, (5) = 3, (6) = 4})

 (3)
> 

Vector([A,B]);

Error, (in Vector) initializer list contains elements of width > 1 and depth > 1

  
> 

Vector(['A','B']);

Vector(2, {(1) = Matrix(2, 2, {(1, 1) = 1, (1, 2) = 2, (2, 1) = 3, (2, 2) = 4}), (2) = Vector(3, {(1) = 2, (2) = 3, (3) = 4})})

 (4)
> 

Matrix(1,2,[[A,B]]);

Error, (in Matrix) this entry is too tall or too short: Vector(3, {(1) = 2, (2) = 3, (3) = 4})

  
> 

Matrix(1,2,[['A','B']])

Matrix(1, 2, {(1, 1) = Matrix(2, 2, {(1, 1) = 1, (1, 2) = 2, (2, 1) = 3, (2, 2) = 4}), (1, 2) = Vector(3, {(1) = 2, (2) = 3, (3) = 4})})

 (5)

Change...

vv 14107
June 03 2017
3 1

J:=int(ln(x)^n,x):
ch := t=ln(x):   IntegrationTools:-Change(J, ch):
simplify(eval(%, ch));

It will work...

vv 14107
June 01 2017
1 0

1. Don't use procedures. Use expressions:

SED2 := 32/3 * ...  

(if you really need procedures, for diff use unapply).

2. To have polynomials extract only the numerators:

SED2 := numer(SED2).

(Maybe you will have to take the derivatives before that, it depends on your intentions).

3. Use ":" instead of ";"  only after everything works OK,

 

Floats?...

vv 14107
June 01 2017
1 0

rand(-2 .. 2)  generates coefficients from the set {-2,-1,0,1,2}.

If you want floats (real numbers), use rand(-2.0 .. 2.0)

Numerically...

vv 14107
May 29 2017
0 1
> 

restart;
> 

Digits:=15;

15

 (1)

 

Solving numerically the recurrence

 

F(0,y) = 0,

F(x,y) = y/(y+1)*F(x-1,y)^((y+1)/y) + 1/(y+1)

 

and approximating  limit( F(n,n),  n = infinity)

 

No chance for a closed form (probably).

 

> 

F:=proc(x,y)
   if x=0 then return 0 fi;
   y/(y+1)*F(x-1,y)^((y+1)/y)+1/(y+1)
end:
> 

F(2,2);

(2/27)*3^(1/2)+1/3

 (2)
> 

F(1000,1000.);

.420589864597856

 (3)
> 

F(10^4,10.^4);

.420522651099579

 (4)
> 

F(10^5,10.^5);

Error, (in F) too many levels of recursion

  
> 

 

 

We must use a non-recursive version:

 

> 

f:=proc(x,y)
local r:=0,a, z:=(y+1)/y, w:=y/(y+1), y1:= 1/(y+1);
to x do r:=  w*r^z+y1 od;
r
end:
> 

f(2,2);

(2/27)*3^(1/2)+1/3

 (5)
> 

seq( [n, f(10^n,10.^n)],n=1..6);

[1, .428422578046699], [2, .421269413312481], [3, .420589864597856], [4, .420522651099579], [5, .420515940874970], [6, .420515270000679]

 (6)
> 

 

Symbolic absolute extrema for two variab...

vv 14107
May 29 2017
2 2


 

> 

restart;

 

The exact (symbolic) absolute min and max.
Two variables. The constraints must define a compact domain.
Based on extrema and solve (which are supposed to work).

> 

absminmax:=proc(f::algebraic, constr::set(relation),candi::name:='_candidates')
local bdr, n,P,i,j,u,s,cand,x,y;
s:=indets([f,constr],name) minus {constants};
if nops(s)<>2 then error "Two indeterminates only" fi;
bdr:=seq(lhs(u)-rhs(u)=0,u=constr):
x,y:=s[];  n:=0:  P:=Vector():
for j to nops(constr) do for i to j-1 do
  s:=solve([bdr[i],bdr[j]],{x,y},explicit);
  for u in [s] do n:=n+1; P(n):=u od;
od od:
for i to nops(constr) do
  extrema(f, {bdr[i]},{x,y}, 's'):
  for u in s do n:=n+1: P(n):=u od;
od:
s:=solve([diff(f,x),diff(f,y)],{x,y},explicit):
for u in [s] do n:=n+1: P(n):=u od:
cand:=select(u -> type(rhs~(u),set(realcons)), {entries(P,nolist)});
cand:=select(u -> is( `and`( eval(constr, u)[]) ), cand);
candi:=cand; if cand={} then return {} fi;
min[index]( [seq(eval(f,u),u=cand)]):  
  s:=`min`=eval(f,cand[%]), cand[%];
max[index]( [seq(eval(f,u),u=cand)]):  
  s, `max`=eval(f,cand[%]), cand[%];
end:

 

Examples

 

 

> 

f1:=2*x^2+4*x*y+2*x-2*y^2+y+4:
> 

constr1:={ y>=(1/21)*x-283/42,
          x<=(3/29)*y+329/58,
          y<=(1/9)*x+131/18,
          x>=-(1/9)*y-113/18 }:
> 

absminmax(f1, constr1);

min = -137657/784, {x = -2267/392, y = 2601/392}, max = 379/2, {x = 13/2, y = 8}

 (1)
> 

_candidates;

{{x = -7, y = 13/2}, {x = 5, y = -13/2}, {x = -2267/392, y = -1745/392}, {x = -2267/392, y = 2601/392}, {x = -11/2, y = -7}, {x = -3/8, y = -1/8}, {x = 13/2, y = 8}}

 (2)
> 

f2:=4*x*y-2*y^2+4*x+y+4:
constr2:={y >= -(1/13)*x-85/13, y <= -(25/14)*x+26/7, y <= (23/12)*x+89/12}:
> 

absminmax(f2, constr2,'ccc');

min = -245, {x = 6, y = -7}, max = 66, {x = -7, y = -6}

 (3)
> 

absminmax(x^2+y^2, {x<=1,y<=3,x>=0,y>=0});

min = 0, {x = 0, y = 0}, max = 10, {x = 1, y = 3}

 (4)
> 

absminmax(x+y, {(x-1)^2+y^2<=1, y>=x^2});

min = 0, {x = 0, y = 0}, max = 2, {x = 1, y = 1}

 (5)
> 

absminmax(exp(x+y), {(x-1)^2+y^2<=1, y>=x^2});

min = 1, {x = 0, y = 0}, max = exp(2), {x = 1, y = 1}

 (6)
> 

absminmax(x+y, {(x-1)^2+y^2<=1, y>=x^2+10});

{}

 (7)
> 

absminmax(x^2-y^2, {(x-1)^2+y^2<=1, y>=x^2});

min = -1/2, {x = 1/2, y = (1/2)*3^(1/2)}, max = 1/4, {x = (1/2)*2^(1/2), y = 1/2}

 (8)
> 

absminmax(x^4-y^2, {(x-1)^2+y^2<=1, y<=x^2});

min = ((1/6)*(54+6*87^(1/2))^(1/3)-1/(54+6*87^(1/2))^(1/3))^4-(1/36)*(-(54+6*87^(1/2))^(4/3)+12*(54+6*87^(1/2))^(2/3)-72*(54+6*87^(1/2))^(1/3)+72*87^(1/2)+612)/(54+6*87^(1/2))^(2/3), {x = (1/6)*(54+6*87^(1/2))^(1/3)-1/(54+6*87^(1/2))^(1/3), y = -(1/6)*((54+6*87^(1/2))^(2/3)*(-(54+6*87^(1/2))^(4/3)+12*(54+6*87^(1/2))^(2/3)-72*(54+6*87^(1/2))^(1/3)+72*87^(1/2)+612))^(1/2)/(54+6*87^(1/2))^(2/3)}, max = 16, {x = 2, y = 0}

 (9)
> 

 


Download absminmax-sent.mw

int...

vv 14107
May 27 2017
1 0

Diff(int(value(z), t), t);

residues...

vv 14107
May 26 2017
1 1


 

> 

restart;
> 

f:=((z^4*exp(2*z)+1)/(z+I)^3) - ((z^3+z)/((z-2*I)*(z-5))) + 8*Pi*exp(1);  # probably

(z^4*exp(2*z)+1)/(z+I)^3-(z^3+z)/((z-2*I)*(z-5))+8*Pi*exp(1)

 (1)
> 

poles:=[-I, 2*I, 5];

[-I, 2*I, 5]

 (2)
> 

C:= z -> is( abs(z-1) < sqrt(11/2) );

proc (z) options operator, arrow; is(abs(z-1) < sqrt(11/2)) end proc

 (3)
> 

2*Pi*I*add( residue(f, z=a), a=select(C, poles) );

(2*I)*Pi*(12/29-(30/29)*I+(102+112*I+(-2-4*I)*((4*I)*(exp(I))^2+5*(exp(I))^2)/(exp(I))^2)/(-(15*I)*(exp(I))^2+3*(exp(I))^2))

 (4)
> 

evalf(%);

4.564531460+58.76512338*I

 (5)
> 

 


 

pdsolve...

vv 14107
May 25 2017
3 0

The fact is that pdsolve is not very competent with initial conditions
and we must do it ad-hoc:

ds:=pdsolve({sys}) ;
           {x(t, s) = _F1(s) sin(t) + _F2(s) cos(t),
             y(t, s) = _F1(s) cos(t) - _F2(s) sin(t)}
ds0:=eval(ds,t=0);
              {x(0, s) = _F2(s), y(0, s) = _F1(s)}
subs(ds0, [ic]);
                 [_F2(s) = a(s), _F1(s) = b(s)]
eval(ds,%);
             {x(t, s) = b(s) sin(t) + a(s) cos(t),
               y(t, s) = b(s) cos(t) - a(s) sin(t)}

 

 

Lists, Vectors , ......

vv 14107
May 25 2017
0 0

For these kinds of maniputations, the vectors, lists, matrices etc were invented.
Example using column vectors:

gamma__w :=10:
z := <2, 1, 0, -1, -4, -7>:
h := <2, 1, 0,  0,  0,  0>:

gamma__w * (h-z);
                   

 

 

 

polynomials...

vv 14107
May 25 2017
1 0

You have a polynomial systems with many parameters.
The result of solve would be huge and in terms of RootOfs.

Just to make an idea, let's eliminate C4 between the first two equations:

poly:=(numer@lhs)~([eq1,eq2,eq3,eq4]):
resultant(poly[1], poly[2], C4);

You will see that even this first step gives a huge polynomial.

For numerical results, give values to parameters and use fsolve.

 

LSSolve...

vv 14107
May 21 2017
1 2

The residuals are much smaller using LSSolve:

eq:=[seq(lhs(EQQ[i]),i=1..36)]:
Digits:=15:
s:=Optimization:-LSSolve(evalf(eq),iterationlimit=100000);

s := [0.960968723053680e-3, [AA[1] = -.347787136052672, AA[2] = 0.510086502851874e-1, AA[3] = 1.16178642934306, AA[4] = .283780970480705, AA[5] = 0.639701227912076e-1, AA[6] = 0.167348852743906e-2, AA[7] = .970111340350232, AA[8] = 0.123075237465221e-1, AA[9] = -2.56229943109674, AA[10] = .119761186424375, AA[11] = 0.514762128059566e-1, AA[12] = 0.653686059489268e-2, AA[13] = 1.17224901537396, AA[14] = 0.123075188025341e-1, AA[15] = -3.10469172766768, AA[16] = .119761196994863, AA[17] = 0.514762189628706e-1, AA[18] = 0.653685044254999e-2, AA[19] = .970010263800571, AA[20] = 0.122243656427024e-1, AA[21] = -2.56235367446998, AA[22] = .119187577004902, AA[23] = 0.512617524055109e-1, AA[24] = 0.636735006533654e-2, AA[25] = 1.17214793882475, AA[26] = 0.122243606766083e-1, AA[27] = -3.10474597103813, AA[28] = .119187587581231, AA[29] = 0.512617585634991e-1, AA[30] = 0.636733990163037e-2, AA[31] = .516890321334600, AA[32] = 0.510086073253494e-1, AA[33] = -1.15838655820066, AA[34] = .283780982288442, AA[35] = 0.639701289588585e-1, AA[36] = 0.167348506238467e-2]]

Digits:=40:
eval( evalf(eq), s[2]);

Unfortunately LSSolve seems to get stuck for Digits>15.

 

 

Because we already know what has to be p...

vv 14107
May 17 2017
1 3

Because we already know what has to be proved, it could be done like this:

 

> 

x1 := arcsin(1/2*(2+(2+(2+2^(1/2))^(1/2))^(1/2))^(1/2))
> 

solve(sin(32*x)=expand(sin(32*x1)),allsolutions);

(1/32)*Pi*_Z1

 (1)
> 

a:=subs( indets(%)[]=n,%);

(1/32)*Pi*n

 (2)
> 

min[index]( [seq(abs(x1 - a), n=1..16)] );

15

 (3)
> 

'x1'=eval(a, n=%);

x1 = (15/32)*Pi

 (4)

 

Looking closer at the bug...

vv 14107
May 17 2017
3 0


 

> 

restart;
> 

sys:= [diff(f(x, y, t, u), u, t)-(diff(f(x, y, t, u), x, y)) = 0,
       diff(f(x, y, t, u), u, u) = 0,
       diff(f(x, y, t, u), u, y) = 0,
       diff(f(x, y, t, u), x, u) = 0,
       diff(f(x, y, t, u), x, x) = 0,
       diff(f(x, y, t, u), y, y, y) = 0];

[diff(diff(f(x, y, t, u), t), u)-(diff(diff(f(x, y, t, u), x), y)) = 0, diff(diff(f(x, y, t, u), u), u) = 0, diff(diff(f(x, y, t, u), u), y) = 0, diff(diff(f(x, y, t, u), u), x) = 0, diff(diff(f(x, y, t, u), x), x) = 0, diff(diff(diff(f(x, y, t, u), y), y), y) = 0]

 (1)
> 

S0:=pdsolve(sys);

{f(x, y, t, u) = (_F3(t)*y+_F4(t))*x+(_F3(t)+_C1)*u+(1/2)*_F7(t)*y^2+_F8(t)*y+_F9(t)}

 (2)
> 

pdetest(S0,sys);  # bug

[diff(_F3(t), t)-_F3(t), 0, 0, 0, 0, 0]

 (3)
> 

S00:=eval(S0, _F3(t)=c*exp(t));

{f(x, y, t, u) = (c*exp(t)*y+_F4(t))*x+(c*exp(t)+_C1)*u+(1/2)*_F7(t)*y^2+_F8(t)*y+_F9(t)}

 (4)
> 

pdetest(S00,sys);  # So, S00 is a solution.  But is it the general one?

[0, 0, 0, 0, 0, 0]

 (5)

Let us try to answer this question.

 

 

> 

#####################################
> 

 

> 

S1:=pdsolve(sys[2..-1]);  # remove the first PDE

{f(x, y, t, u) = (1/2)*(2*_F5(t)+y*(_F3(t)*y+2*_F4(t)))*x+_F8(t)*u+(1/2)*_F9(t)*y^2+_F10(t)*y+_F11(t)}

 (6)
> 

pdetest(S1,sys[2..-1]); #ckeck OK

[0, 0, 0, 0, 0]

 (7)
> 

pdetest(S1,sys);

[diff(_F8(t), t)-_F3(t)*y-_F4(t), 0, 0, 0, 0, 0]

 (8)
> 

solve(%,_F4(t));

{_F4(t) = -_F3(t)*y+diff(_F8(t), t)}

 (9)
> 

S2:=eval(S1,%);

{f(x, y, t, u) = (1/2)*(2*_F5(t)+y*(-_F3(t)*y+2*(diff(_F8(t), t))))*x+_F8(t)*u+(1/2)*_F9(t)*y^2+_F10(t)*y+_F11(t)}

 (10)
> 

# S2 should be the general solution of the initial system
> 

pdetest(S2,sys);

[_F3(t)*y, 0, 0, 0, 0, 0]

 (11)
> 

# But it is not. Conclusion: the solution S1 of the system S[2..-1] is not the general one.
# So, probably S00 is not general too.
# It should be simple to check this by hand.
> 

###########################################
> 

 


 

Download pde-bug.mw

Edit.  I made a mistake when solving for _F4.
_F3  must be necessarily 0 (due to the presence of y)
So, actually S2 must be

So, indeed, S00 is not general (for example f = t*u+x*y is not contained in S00), but S1 could be after all correct; it would be interesting to check this by hand.
 

 

 

proc...

vv 14107
May 15 2017
2 1

 

subs[eval](F = ( (x,y)->f(x)+g(y) ), F(x,x) );

(Of course [eval] is not mandatory, just to obtain th result evaluated.)

First 89 90 91 92 93 94 95 Last Page 91 of 121