vv

1734 Reputation

9 Badges

1 years, 121 days

MaplePrimes Activity


These are answers submitted by vv


 

The 3 circles do not have a common point of intersection.
But one point in the intersection of (eq1 and eq2)  
is very close to one poin in the intersection of (eq2 and eq3).
You can find this point e.g. minimizing eq1^2+eq2^2+eq3^2

 

restart;

eq1:= x^2+y^2-6600*x-4400*y+15717900;
eq2:= x^2+y^2-6820*x-4180*y+15989800;
eq3:= 1.966975428*10^7+x^2+y^2-7480*x-4840*y;

x^2+y^2-6600*x-4400*y+15717900

 

x^2+y^2-6820*x-4180*y+15989800

 

19669754.28+x^2+y^2-7480*x-4840*y

(1)

s12:={solve([eq1,eq2],explicit)};evalf(%);

{{x = 74095/22-(5/22)*27727^(1/2), y = 46905/22-(5/22)*27727^(1/2)}, {x = 74095/22+(5/22)*27727^(1/2), y = 46905/22+(5/22)*27727^(1/2)}}

 

{{x = 3330.110394, y = 2094.201304}, {x = 3405.798696, y = 2169.889606}}

(2)

s23:={solve([eq2,eq3],explicit)};

{{x = 3405.798698, y = 2169.889605}, {x = 3489.889605, y = 2085.798698}}

(3)

s12 intersect s23;

{}

(4)

minimize( eq1^2+eq2^2+eq3^2, location );

-0.2084e19, {[{x = 3405.798701, y = 2169.889609}, -0.2084e19]}

(5)

 

EDIT. The value of the minimum given by minimize is of course wrong (bug!).

loc:=%[2][][1];
               {x = 3405.798701, y = 2169.889609}
eval([eq1,eq2,eq3],loc);
                       [0.01, 0., -0.01]

 

 

EQ:=convert(lhs(eq1),set) union convert(lhs(eq2),set):
sol:=solve(EQ):
eval([fn1,gn1],sol);

N.B. You should avoid evalm and use "+", ".".

Use:

[A$2] .~ LM;

or

(A$2) .~ LM;

 

 

GetCoeff:=proc(f,u)
  local x:=indets(f,name);
  if u=0 or u=1 then return eval(f,x=~0) fi; 
  eval(foldl(coeff,f,op(u)), x=~0);
end proc:

 

It is possible to represent (and "see")  only surfaces. You cannot see in the interior of a 3-dimensional set in R^3 (between any interior point and your eye there exist an infinity of other interior points obfuscating the first one).

Today you have problems easier by hand than with Maple :-)

f = (x+y)/(x+y+z), g = (x+z)/(x+y+z), h = (y+z)/(x+y+z)  have the same integral (by symmetry).

==> int(f) = 1/3*int(f+g+h) = 1/3*int(2) = 1/3*2*(4*Pi/3)/8 = Pi/9.

Edit.
A real challenge for any CAS would be to compute the integral of
(x^2017+y^2017)/(x^2017+y^2017+z^2017)
in the first octant of the unit ball. [=Pi/9]

dxdy  :=   y=-sqrt(1-x^2) .. sqrt(1-x^2), x=-1..1 ;


int(sin(x^2) *cos(x^2) , dxdy) =

1/2*int( sin(x^2+y^2) + sin(x^2 - y^2), dxdy ) =

1/2*int( sin(x^2+y^2) , dxdy)  + 0 =

1/2* int( sin(r^2)*r, r=0..1, t=0..2*Pi ) =

1/2*(1 - cos(1))*Pi =

Pi * sin(1/2)^2

 

It's a very difficult task to compute u3 even at a single point.
For example, u3(0) can be simplified to the triple integral of
f := -(2.533317577*10^7*I)*z*exp((-.64-12.56637062*I)*x^2-(7.853982892*10^7*I)*z^2-(25.13274124*I)*y^2)*y*x*BesselJ(0, 25.13274123*y*x)*BesselJ(0, 62831.85308*z*y)

But f is very oscillatory and to compute its triple integral over [0,1]^3 is almost impossible even using MonteCarlo method.
Probably some special methods will be needed.
Just look at this plot:
plot( eval(Re(f),[x=1/2,y=1/2]),z=0..1 );


 

P.S. Was this integral invented, or it appears in a concrete problem?

eq:=proc(L::listlist,x::name,y::name,z::name)
  primpart(LinearAlgebra:-Determinant(<Matrix([[x,y,z],L[]])|<1,1,1,1>>))=0
end:

map(eq,L,x,y,z);

 

TriOnUnitSphere:=proc(A::Vector[column](3),B::Vector[column](3),C::Vector[column](3)) #A,B,C = directions
uses LinearAlgebra,plots,plottools;
local a:=Normalize(A,2),b:=Normalize(B,2),c:=Normalize(C,2),u,v;
display(sphere(),  plot3d(1.01*Normalize(c+v*(a+u*(b-a)-c),2), u=0..1,v=0..1,color=red,style=surface,_rest),
pointplot3d([1.01*a,1.01*b,1.01*c],symbolsize=20,color=blue,symbol=solidsphere)  )
end:

TriOnUnitSphere(<1,2,3>, <-1,1,-2>, <-1,-1,3>);

TriOnUnitSphere(<1,0,0>, <-10,1,0>, <0,1,5>, grid=[200,200]);

 

with(LinearAlgebra):with(plots):with(plottools):
A1:=<0,0,1>: A2:=<0,3/5,4/5>: A3:=<4/5,3/5,0>: A4:=<5/13,0,12/13>:
seg:= (a,b)->spacecurve(Normalize(a+t*(b-a),2),t=0..1, thickness=3, color=red):
display(sphere([0,0,0],1),seg(A1,A2),seg(A2,A3),seg(A3,A4),seg(A1,A4));

 


 

# workaround

eq:=(2*y(x)*diff(diff(y(x),x),x)-diff(y(x),x)^2)^3+32*diff(diff(y(x),x),x)*(x*diff(diff(y(x),x),x)-diff(y(x),x))^3 = 0;

(2*y(x)*(diff(diff(y(x), x), x))-(diff(y(x), x))^2)^3+32*(diff(diff(y(x), x), x))*(x*(diff(diff(y(x), x), x))-(diff(y(x), x)))^3 = 0

(1)

solve(eq,  diff(diff(y(x),x),x));

RootOf(32*x^3*_Z^4+(-96*x^2*(diff(y(x), x))+8*y(x)^3)*_Z^3+(-12*(diff(y(x), x))^2*y(x)^2+96*x*(diff(y(x), x))^2)*_Z^2+(6*(diff(y(x), x))^4*y(x)-32*(diff(y(x), x))^3)*_Z-(diff(y(x), x))^6)

(2)

eq1:=diff(diff(y(x),x),x)=%

diff(diff(y(x), x), x) = RootOf(32*x^3*_Z^4+(-96*x^2*(diff(y(x), x))+8*y(x)^3)*_Z^3+(-12*(diff(y(x), x))^2*y(x)^2+96*x*(diff(y(x), x))^2)*_Z^2+(6*(diff(y(x), x))^4*y(x)-32*(diff(y(x), x))^3)*_Z-(diff(y(x), x))^6)

(3)

sol:=dsolve( eq1);

y(x) = (2/3)*_C1^2*x^2*3^(2/3)/(_C2*(_C1*(3^(1/2)*(_C1*(27*_C1*_C2-64)/_C2)^(1/2)+9*_C1)*_C2^2)^(1/3))+(1/6)*_C1*(_C1*(3^(1/2)*(_C1*(27*_C1*_C2-64)/_C2)^(1/2)+9*_C1)*_C2^2)^(1/3)*x^2*3^(1/3)/_C2^2+(1/4)*_C1^2*x^2/_C2+_C1*x+_C2

(4)

odetest(sol,eq);

0

(5)

DEtools:-odeadvisor(eq1);

[NONE]

(6)

You can use e.g.

binaryvariables = indets(L, specindex(integer,x));

where L is a list containing the constraints (and the function).
Or, you may generate their sequence using seq.
 

f:=(x,y) ->exp(-x^2-y^2):
L:=4: n:=64:
r1:=rand(0..1.0): r:=rand(-1..1.0):
a:=Vector(n,L*r): b:=Vector(n,L*r):
c:=Vector(n,0.3*r):
da:=Vector(n, i->0.5*''r1()''):
P:=proc(aa)
  global n,a,b,c,da,h,L;
  a:=a+eval~(da);
  map[inplace](frem,a,L);
  h:=add(c[i]*f(x+a[i], y+b[i]),i=1..n):
  plot3d(h, x=-L..L,y=-L..L, color="Aquamarine", scaling=constrained, view=-2..2, axes=none, style=surface, lightmodel=light1);
end:
plots:-display(seq(P(aa),aa=1..50),insequence=true);
## or ...
Explore(P(aa), aa=0..1, animate,loop,size=[800,800]);

 

There are many possible recurrences. The second one is yours.

 

restart;

u:=n->(1-x^2)^n/n!;

proc (n) options operator, arrow; (-x^2+1)^n/factorial(n) end proc

(1)

rec:=diff(f[n](x),x,x) + a*f[n-1](x) + b*x^2* f[n-2](x);

diff(diff(f[n](x), x), x)+a*f[n-1](x)+b*x^2*f[n-2](x)

(2)

diff(u(n),x,x) + a*u(n-1) + b*x^2* u(n-2): simplify(%/u(n-2));

(((b+4)*n-a-b-2)*x^2+a-2)/(n-1)

(3)

coeffs(expand(%),x):

solve({%},{a,b});

{a = 2, b = -4}

(4)

eval(rec,%)=0;

diff(diff(f[n](x), x), x)+2*f[n-1](x)-4*x^2*f[n-2](x) = 0

(5)

############################

rec:=diff(f[n](x),x,x) + a*f[n-1](x) + b* f[n-2](x);

diff(diff(f[n](x), x), x)+a*f[n-1](x)+b*f[n-2](x)

(6)

diff(u(n),x,x) + a*u(n-1) + b* u(n-2): simplify(%/u(n-2));

((4*x^2+b)*n+(-a-2)*x^2+a-b-2)/(n-1)

(7)

coeffs(expand(%),x):

solve({%},{a,b});

{a = 4*n-2, b = -4}

(8)

eval(rec,%)=0;

diff(diff(f[n](x), x), x)+(4*n-2)*f[n-1](x)-4*f[n-2](x) = 0

(9)

#######################

 


 

Download req-int.mw

1 2 3 4 5 6 7 Last Page 1 of 24