Kitonum

13999 Reputation

24 Badges

11 years, 189 days

MaplePrimes Activity


These are answers submitted by Kitonum

The logic is very simple. The   plots:-implicitplot  command works by calculating the values of a function of two variables  F(x,y)  on a grid. For example, if  F(x1,y1)<0  and  F(x2,y2)>0 , then on the segment connecting these points there exists a point such that  F(x0,y0)=0  (of course we assume that the function  F  is continuous.
Your function  F(r,theta)=r-cos(theta)  in polar coordinates . It is equal to 0 at points of the circle with center  (x,y)=(1/2, 0)  and radius 1/2. You are plotting an region outside this circle. I replaced  r=1/3  in your code with  r=1.5  so that the whole circle is visible. For unknown reason, Maple does not paint over the region to the left of the Oy axis. For the correct plotting, we can use the Cartesian coordinates:


 

restart;

plots:-implicitplot(r>=cos(theta), r = 0 .. 1.5, theta = 0 .. 2*Pi, filledregions, coords = polar, numpoints=5000,  axiscoordinates = polar);

plots:-implicitplot((x-1/2)^2+y^2>=1/4, x=-0.5..1.5,y=-1..1, filledregions, numpoints=5000, scaling=constrained);

 

 

 


Addition. If there is no sign change, then the  plots:-implicitplot  command fails. Compare (obviously abs(x-y)=0 is equivalent to x-y=0) .

plots:-implicitplot(abs(x-y)=0, x=-1..1,y=-1..1);
plots:-implicitplot(x-y=0, x=-1..1,y=-1..1);

Download region_polar.mw

 

An easy way is to write both equations in one system. Here is an example:

 

restart;         
Sol:=dsolve({diff(P(x),x)=2*x,diff(V(x),x)=cos(x),P(0)=0,V(0)=0},numeric):
plots:-odeplot(Sol,[P(x),V(x)], x=0..10);
 

 

 


Edit. If you already have separate solutions for these equations (the same equations), then you can do so

restart; 	
Sol1:=dsolve({diff(P(x),x)=2*x,P(0)=0},numeric):
Sol2:=dsolve({diff(V(x),x)=cos(x),V(0)=0},numeric):
plot([t->eval(P(x),Sol1(t)), t->eval(V(x),Sol2(t)),0..10]);

The result is the same as above.

 

Download PV.mw

We can make the animation much smoother if we provide linear interpolation between successive frames:

restart;     
SOLNSuy[1, 1] := 2.5872902469406659197*10^(-20)-.65694549571241255901*y
+1.9708364871372376767*y^2-1.3138909914248251176*y^3-1.6010739356637904911*10^(-19)*y^4:     
SOLNSuy[2, 1] := -4.002204462000*10^(-20)-1.7879176897079605225*y
+5.3637530691192141414*y^2-3.5758353794044226250*y^3-6.8309939211286845440*10^(-12)*y^4:     
SOLNSuy[3, 1] := -1.1953264450000*10^(-19)-3.2481690589079594122*y
+9.7445071767154794599*y^2-6.4963381177952273213*y^3-1.2292726248071398400*10^(-11)*y^4:    
SOLNSuy[4, 1] := -2.6720465500000*10^(-19)-4.9239979672954025921*y
+14.771993901873204315*y^2-9.8479959345587718955*y^3-1.9029826928878336000*10^(-11)*y^4:     
SOLNSuy[5, 1] := 3.416928541000*10^(-20)-6.7268498492441931137*y
+20.180549547714413714*y^2-13.453699698443639810*y^3-2.6580790570532587008*10^(-11)*y^4:     
SOLNSuy[6, 1] := -2.554122292000*10^(-20)-8.5884528335125514887*y
+25.765358500514014457*y^2-17.176905666966875698*y^3-3.4587270427710613504*10^(-11)*y^4:     
SOLNSuy[7, 1] := -9.206107680000*10^(-20)-10.456823708331499352*y
+31.370471124965259849*y^2-20.913647416590986491*y^3-4.2774005353527132160*10^(-11)*y^4:     
SOLNSuy[8, 1] := 1.9644186790000*10^(-19)-12.293003938471349390*y
+36.879011815379230436*y^2-24.586007876856948223*y^3-5.0932823222176363520*10^(-11)*y^4:    
SOLNSuy[9, 1] := -3.775112769000*10^(-19)-14.068404975282556550*y
+42.205214925807397100*y^2-28.136809950465931724*y^3-5.8908824448577377280*10^(-11)*y^4:     
SOLNSuy[10, 1] := 1.146281780000*10^(-19)-15.762658869974768890*y
+47.287976609878780960*y^2-31.525317739837422477*y^3-6.6589592851037286400*10^(-11)*y^4:

t:=a-floor(a):
F:=a->(1-t)*SOLNSuy[floor(a), 1]+t*SOLNSuy[floor(a)+1, 1]:
     
plots[animate](plot, [F(a), y=0..1], a = 1 .. 10, frames=50, size=[800,400]);

            

If it is necessary to introduce a new function (let's say with the name  g )  such that  g(zeta)=f(x)  under the condition  zeta=x/a , then it is obvious that  g:=zeta->f(a*zeta)
So the new function  g  will be

f:=x->A1*sin(k*x)+A2*cos(k*x)+A3*sinh(k*x)+A4*cosh(k*x):
g:=unapply(f(a*zeta), zeta);

  

 

If we need to select terms of a certain degree (relative to some variable) in a polynomial, then this can be done as in this example:

P:= y^2 - x^2/y:
n:=degree(P,x);
select(t->degree(t,x)=n, P);

                        

Use the Physics package for this (we get the result in dot-notation):

restart;
T__ln := (1/2)*m__ws*v__l+(1/2)*I__ws*(diff(varphi__l(t), t))^2+((1/2*(m__zs+m__zp+m__wp))*v__l+(1/2*(I__zs+I__zp+I__wp))*(diff(varphi__l(t), t))^2)*z__11*eta/z__12+((1/2*(m__wk+m__zpk+m__k+m__zk))*v__l+(1/2*(I__wk+I__zpk+I__k+I__zk))*(diff(varphi__l(t), t))^2)*z__21*eta/z__22;
with(Physics):
diff(T__ln, diff(varphi__l(t), t));

               

I advise you to master the basic Maple commands, which include the  plot  command for plotting of functions. The  discont  option allows you to avoid vertical segments at the break points of the function. Maple automatically displays the values of this function at break points in the form of small squares:

plot(floor(x)+3, x=-4..3, color=red, discont, scaling=constrained, size=[500,500]);

                    

 

Note that if  is numeric, x mod 3 returns a result only for an integer  x . Otherwise, an error occurs. Therefore, your plot will consist of separate points corresponding to integer values of  x . There is a special syntax for plotting individual points (see below). You can present these results in a matrix or table (see below).


 

plot([seq([x,x mod 3],x=-5..5)],color=red, style=point, symbol=solidcircle,scaling=constrained,size=[1000,300]);
interface(rtablesize=100):
A:=<`<|>`(<x,"x mod 3">,seq(<x,x mod 3>,x=-5..5))>;
DocumentTools:-Tabulate(A, width=100):

 

Matrix(2, 12, {(1, 1) = x, (1, 2) = -5, (1, 3) = -4, (1, 4) = -3, (1, 5) = -2, (1, 6) = -1, (1, 7) = 0, (1, 8) = 1, (1, 9) = 2, (1, 10) = 3, (1, 11) = 4, (1, 12) = 5, (2, 1) = "x mod 3", (2, 2) = 1, (2, 3) = 2, (2, 4) = 0, (2, 5) = 1, (2, 6) = 2, (2, 7) = 0, (2, 8) = 1, (2, 9) = 2, (2, 10) = 0, (2, 11) = 1, (2, 12) = 2})

(1)


Download mod.mw

In the figure, the red broken  AEC  is the above border of the first assembly, and the red broken  ABC  is the above border of the second one. The blue line  AC  is a straight line. Using the  geometry  package, we show that the area of the quadrangle  ABCE  is exactly  1 .

restart:
with(plottools): with(plots): with(geometry):  
P1:=polygon([[0,0],[8,0],[8,3]], color="Pink", style=surface):
P2:=polygon([[8,0],[8,1],[10,1],[10,2],[13,2],[13,0]], color=green, style=surface):
P3:=polygon([[8,1],[8,3],[13,3],[13,2],[10,2],[10,1]], color="Yellow", style=surface):
P4:=polygon([[8,3],[13,5],[13,3]], color="Green", style=surface):
L:=plottools:-line([0,0],[13,5], color=blue):
L1:=curve([[0,0],[8,3],[13,5]],color=red,thickness=2):
L2:=curve([[0,0],[5,2],[13,5]],color=red,thickness=2):
A:=[0,0]: B:=[5,2]: C:=[13,5]: E:=[8,3]:
T:=textplot([[A[],"A"],[B[],"B"],[C[],"C"],[E[],"E"]],font=[times,20],align={above,left}):
P:=pointplot([A,B,C,E], symbol=solidcircle, color=red, symbolsize=8):
display(L,L1,L2, P1,P2,P3,P4,T,P, scaling=constrained, size=[1100,450], axis[1] = [gridlines = [13,color=grey,thickness=0]], axis[2] = [gridlines = [5,color=grey,thickness=0]], transparency=0.5);
area(triangle(t1,[point(AA,A),point(BB,B),point(EE,E)]))+area(triangle(t2,[point(CC,C),point(BB,B),point(EE,E)]));

        

I agree that this is a serious bug. It seems that in the second case Maple first rounds each of the operands to 2 digits, and only then performs the calculation. Therefore, we get  1000-1000 = 0 .
Here is a workaround:
restart:  
P1 := 1007:
P3 := 1014.1:
P3 - P1:
evalf(%, 2);
                                
7.1

A useful conclusion: first you need to do all the calculations with maximum accuracy (preferably symbolically), and only at the end, round to the desired number of digits.

In your condition

, instead of  x , there should be a constant. I took  x=b :


 

"restart:Digits :=15: upsilon:=0.3:E(x):=E0*((x)/((b)))^(beta):rho(x):=rho0*((x)/((b)))^(beta):alpha(x):=alpha0*((x)/((b)))^(beta):a:=0.2:b:=1:omega:=100:E0:=390e9:rho0:=3900:T(x):=Ta+(Tb-Ta)/(ln(b/(a)))*(ln(x)-ln(a)):Ta:=373:Tb:=273:upsilon:=0.25:alpha0:=7e-6:  h(x):=(1-n*(x/(b)))^(k):n:=0.415196:k:=3:beta:=1:    dsys5 := {(1/(b))*( diff(u(x),x,x) )+(1/(b*h(x))*(diff(h(x),x))+1/(b*E(x))*(diff(E(x),x))+1/(b*(x)))*(diff(u(x),x))+((upsilon)/((b^(2)*x))*1/(h(x))*(diff(h(x),x))-1/((b*x)^(2))+(upsilon)/(b^(2)*(x))*1/(E(x))*(diff(E(x),x)))*b*u(x)+(1+upsilon)*((rho(x)*x*b*(omega^(2)))/(E(x))*(1-upsilon)-(alpha(x)*Ta)/(b)*(diff(T(x),x))-((diff(alpha(x),x))/(b)+(alpha(x)*diff(E(x),x))/(b*E(x))+(alpha(x)*diff(h(x),x))/(b*h(x)))*Ta*T(x) ),u(a) = 0,(E(b))/((1-upsilon^(2)))*(D^((1))(u)(b)+upsilon/(b)*D^((0))(u)(b))-(E(b)*alpha(b)*T(b)*Ta)/((1-upsilon^())) =-1}:dsol5 := dsolve(dsys5,abserr=1e-1, 'maxmesh'=900, numeric, method=bvp[middefer],output=listprocedure):fy := eval(u(x),dsol5)"

proc (x) local res, data, solnproc, `u(x)`, outpoint; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](x) else outpoint := evalf(x) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(8, {(1) = .2, (2) = .2973512592166774, (3) = .4046598589973246, (4) = .5253877794811062, (5) = .6502960591623015, (6) = .7736319028878446, (7) = .8952969100357793, (8) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(8, 2, {(1, 1) = .0, (1, 2) = .733488879983426, (2, 1) = 0.6165197085959364e-1, (2, 2) = .5330991714843657, (3, 1) = .1180091558271178, (3, 2) = .5172767551930041, (4, 1) = .18269975228050025, (4, 2) = .554399062532572, (5, 1) = .255347334921374, (5, 2) = .6088157832133286, (6, 1) = .33396482109823816, (6, 2) = .6660364219868107, (7, 1) = .4183769220601095, (7, 2) = .7215786851118666, (8, 1) = .4963025525973446, (8, 2) = .76692811184826}, datatype = float[8], order = C_order); YP := Matrix(8, 2, {(1, 1) = .733488879983426, (1, 2) = -4.437310575901379, (2, 1) = .5330991714843657, (2, 2) = -.6408082738085314, (3, 1) = .5172767551930041, (3, 2) = .16264757217660253, (4, 1) = .554399062532572, (4, 2) = .3928795175533184, (5, 1) = .6088157832133286, (5, 2) = .45496331913237875, (6, 1) = .6660364219868107, (6, 2) = .461339583360867, (7, 1) = .7215786851118666, (7, 2) = .4443915660014327, (8, 1) = .76692811184826, (8, 2) = .41749625147947433}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(8, {(1) = .2, (2) = .2973512592166774, (3) = .4046598589973246, (4) = .5253877794811062, (5) = .6502960591623015, (6) = .7736319028878446, (7) = .8952969100357793, (8) = 1.0}, datatype = float[8], order = C_order); Y := Matrix(8, 2, {(1, 1) = .0, (1, 2) = 0.9052597788842265e-3, (2, 1) = -0.17533040428448915e-2, (2, 2) = 0.11363665769424824e-1, (3, 1) = -0.19873923968831836e-2, (3, 2) = 0.7223177987320639e-2, (4, 1) = -0.17529391294966445e-2, (4, 2) = 0.3957014218449073e-2, (5, 1) = -0.14983216331688753e-2, (5, 2) = 0.21532793401388695e-2, (6, 1) = -0.131401163357453e-2, (6, 2) = 0.11818256841032992e-2, (7, 1) = -0.11878674118888655e-2, (7, 2) = 0.5977854547565276e-3, (8, 1) = -0.1108948738838933e-2, (8, 2) = 0.2772371847096222e-3}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[8] elif outpoint = "order" then return 2 elif outpoint = "error" then return 0.113636657694248244e-1 elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [2, 8, [u(x), diff(u(x), x)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[8] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[8] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(2, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(8, 2, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(2, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(8, 2, X, Y, outpoint, yout, L, V) end if; [x = outpoint, seq('[u(x), diff(u(x), x)]'[i] = yout[i], i = 1 .. 2)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[8] elif outpoint = "order" then return 2 elif outpoint = "error" then return 0.113636657694248244e-1 elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [2, 8, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[8] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[8] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(2, {(1) = .0, (2) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(8, 2, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(2, {(1) = 0., (2) = 0.}); `dsolve/numeric/hermite`(8, 2, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 2)] end proc, (2) = Array(1..3, {(1) = 18446745813876114814, (2) = 18446745813876115254, (3) = 18446745813876115430}), (3) = [x, u(x), diff(u(x), x)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(x) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(x) else `u(x)` := pointto(data[2][2]); return ('`u(x)`')(x) end if end if; try res := solnproc(outpoint); res[2] catch: error  end try end proc

(1)

plots:-odeplot(fy, [x, u(x)], x = a .. b, color = red)

 

``


 

Download 2023_new.mw

restart;
pointset := [A, [0,0], B, [0,v_B], C, [u_C, v_C], DD, [x_DD, y_DD], OO, [x_OO, y_OO]];
for i from 1 to nops(pointset) by 2 do
assign(pointset[i],pointset[i+1]);
od:

 A, B, C, DD, OO;  
                [0, 0], [0, v_B], [u_C, v_C], [x_DD, y_DD], [x_OO, y_OO]


Take a look at what your code does:

restart;
pointset := [A, [0,0], B, [0,v_B], C, [u_C, v_C], DD, [x_DD, y_DD], OO, [x_OO, y_OO]]:
for i from 1 to nops(pointset) by 2 do
pointset[i]:=pointset[i+1];
od:
A, B, C, DD, OO;
seq(pointset[i], i=[1,3,5,7,9]);

                                          A, B, C, DD, OO
                    [0, 0], [0, v_B], [u_C, v_C], [x_DD, y_DD], [x_OO, y_OO]

The difference between  X:=Y and assign(X,Y)  is that in the first case  X  is not evaluated.
I replaced the names  D  and  O  with  DD  and  OO  because  D  and  O  are protected.

Of course, it’s most convenient to use a multiple assignment in your problem:

A, B, C, DD, OO :=   [0, 0], [0, v_B], [u_C, v_C], [x_DD, y_DD], [x_OO, y_OO];

restart;
#FIRST STEP: Lets find the recursive formula for 
d:=proc(m)
option remember;
if m=1 then return 1 else
expand(-add(xi[i]*d(m-i+1),i=2..m)) fi;
end proc:
#lets check
for i from 1 to 4 do
d(i);
end do; 


The code for the matrix:

DD:=m->Matrix(m,m, (i,j)->`if`(i>j,0,d[j-i+1])):
DD(4);

 

Just solve these equations together by taking different names for the parameters. If the number of solutions is infinite, then the lines coincide. If Maple returns a unique solution, then the lines intersect. If Maple returns NULL (no solutions), then the lines are parallel:

restart;
L1:=<2-4*t,5+6*t>: L2:=<-6-12*s,17+18*s>: L3:=<20*t,-30*t>:
solve({L1[1]=L2[1],L1[2]=L2[2]}); # The lines are identical
solve({L2[1]=L3[1],L2[2]=L3[2]}); # The lines are parallel

  Output:                                   {s = s, t = 2+3*s}
                                                         NULL


Additionally, in the first example, a relationship is established between the parameters. In 3D, the solution is about the same. The straight lines will be skew if there are no solutions and the direction vectors (coefficients in front of the parameters) are not collinear. Writing a procedure for all of these is trivial.                           

You have a linear system with unknowns  [A,B,C,E,L,F], so you can solve it using  LinearAlgebra:-LinearSolve  command. Maple does not return the whole solution (due to its bulkiness), but you can get the values of individual unknowns (in the file  below the first unknown is returned):

 

Download sys.mw

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