May be as follows:

solve({x=z,y=z, z=t}, {x,y,z});  # t is the parameter

             {x = t, y = t, z = t}

If not, please give the full text of the problem.

To be more precise, the planets revolve not in circles, but in ellipses, and the Sun is at one focus. So first we take the polar equation of the ellipse (the Sun is in the pole), and then all the constructions using parametric equations, based on the original polar equations. 

For example, all parameters are arbitrary and have no relation to the actual planets. The plane of rotation of the red planet taken at an angle of 18 degrees with respect to the plane of the blue planet. Rotation period of the blue planet taken 2 times more than the red planet. Initial phases are also different.

As a result of the code, you can see two animations. The first animation shows the motion of the two planets in their orbits, the second one also shows the planes of the orbits:

restart;
r1 := (3*(1-0.4^2))/(1-0.4*cos(t)):
r2 := (2*(1-0.5^2))/(1-0.5*cos(2*t-(1/4)*Pi)):
R := <cos(phi), 0, -sin(phi); 0, 1, 0; sin(phi), 0, cos(phi)>:
phi := (1/10)*Pi:
L1 := [r1*cos(t), r1*sin(t), 0]:
L2 := convert(R.`<,>`(r2*cos(2*t-(1/4)*Pi), r2*sin(2*t-(1/4)*Pi), 0), list): T1 := plots[spacecurve](L1, t = 0 .. 2*Pi, linestyle = 3, color = blue, thickness = 2):
T2 := plots[spacecurve](L2, t = 0 .. 2*Pi, linestyle = 3, color = red, thickness = 2):
Sun := plottools[sphere]([0, 0, 0], 0.2, style = surface, color = yellow):
S1 := plot3d([u*cos(t), u*sin(t), 0], u = 0 .. r1, t = 0 .. 2*Pi, style = surface, color = "LightBlue", transparency = 0.5):
S2 := plot3d(convert(R.<v*cos(2*t-(1/4)*Pi), v*sin(2*t-(1/4)*Pi), 0>, list), v = 0 .. r2, t = 0 .. 2*Pi, style = surface, color = pink):
P1 := proc (x, y) plots[display](plottools[sphere](x, y, style = surface, color = blue)) end proc:
P2 := proc (x, y) plots[display](plottools[sphere](x, y, style = surface, color = red)) end proc:
Planet1 := plots[animate](P1, [L1, 0.15], t = 0 .. 2*Pi, frames = 50):
Planet2 := plots[animate](P2, [L2, 0.1], t = 0 .. 2*Pi, frames = 50):

plots[display](T1, T2, Sun, Planet1, Planet2, scaling = constrained, axes = framed, orientation = [65, 75], lightmodel = light1);
plots[display](T1, T2, S1, S2, Sun, Planet1, Planet2, scaling = constrained, axes = framed, orientation = [65, 75], lightmodel = light1);

evalf(map(log10, [25,5,1,10,4,20]));

map(t->10^t, %);

[1.397940008, .6989700041, 0., 1., .6020599914, 1.301029996]

[24.99999996, 4.999999997, 1., 10., 4.000000001, 20.00000002]

Maybe you mean the plane domain  sqrt(x)<y<=1/sqrt(2) ?

Plotting in Maple 12 classic:

restart;

A:=plot(sqrt(x), x=-0.2..0.7, color=blue, thickness=2):

B:=plot(1/sqrt(2), x=-0.2..0.7, color=red, thickness=2):

A1:=plot(sqrt(x), x=0..1/2, color=white, filled=true):

B1:=plot(1/sqrt(2), x=0..1/2, color=yellow, filled=true):

B2:=plot(1/sqrt(2), x=0..1/2, color=yellow, thickness=3):

C:=plot([1/2, t, t=0..1/sqrt(2)], linestyle=3, color=black):

C1:=plots[textplot]([1/2,-0.02,1/2], font=[TIMES,ROMAN,16]):

plots[display](B2, A, B, A1, B1, C, C1, scaling=constrained);

Of course, manual axes are very laborious for a single application. But it is easy to write a procedure in which to provide any options. 

restart;

export_plot_options := axes = none, size = [850, 850]:
points := [seq([seq(exp(-(x^2+y^2)*(1/100)), x = -10 .. 10.0)], y = -10 .. 10.0)]:

A := plots:-listdensityplot(points, export_plot_options):
B := plots[arrow]([.43-2, .43], [27, 0], width = 0.4e-1, head_length = [0.4e-1, relative], head_width = [15, relative]):
C := plots[arrow]([.43, .43-2], [0, 27], width = 0.4e-1, head_length = [0.4e-1, relative], head_width = [15, relative]):
Ticks := seq(plottools[line]([2*i, .43], [2*i, .43-.3], thickness = 3), i = 0 .. 11), seq(plottools[line]([.43-.3, 2*i], [.43, 2*i], thickness = 3), i = 0 .. 11):
SubTicks := seq(plottools[line]([i, .43], [i, .43-.1], thickness = 3), i = -1 .. 23), seq(plottools[line]([.43-.1, i], [.43, i], thickness = 3), i = -1 .. 23):
Text := plots[textplot]([seq([2*i, .43-.8, 2*i], i = 1 .. 11), [25, .43-.8, x], seq([.43-1, 2*i, 2*i], i = 1 .. 11), [.43-1, 25, y]], font = [TIMES, roman, 30]):

plots[display](A, B, C, Ticks, SubTicks, Text);

I fixed a few mistakes in your code.. Symbolic solution is very cumbersome, so I show only the approximate one and its plot:

restart;

eq:=a11*diff(phi(x),x,x,x,x)+a22*diff(phi(x),x,x)+a33*phi(x):

inc:=phi(a)=sigma1, phi(-a)=sigma1, D(phi)(a)=0, D(phi)(-a)=0:

a11:=2.731e-10:

a22:=-1.651e-9:

a33:=3.09027e-10:

a:=35.714:

sigma1:=200e6:

dsolve({eq, inc});

evalf(%);

assign(%);

plot(phi(x), x=-a-1..a+1, thickness=2);

Detailed explanations about this methods see  here  and  here. 

Solutions your task  (y=1/x^(1/4)  implies  x=1/y^4):

Pi*int((1/x^(1/4))^2-1^2, x=1/16..1);  # Disk integration (integrate along  x-axis)

2*Pi*int(y*(1/y^4)-y*(1/16), y=1..2);   # Shell integration (integrate along  y-axis)

                                  9/16*Pi

                                  9/16*Pi

Vizualization your region (Maple 12 classic):

A:=plot([1/x^(1/4), [1/16,t,t=-1/4..3], 1], x=-1/4..3/2, -1/4..3, thickness=2, color=black):

B:=plot(1/x^(1/4), x=1/16..1, color=yellow, filled):

C:=plot(1, x=1/16..1, color=white, filled):

plots[display](A, C, B);

It is simpler to find the coordinates of all of displaced vertices. For clarity, I have increased the all displacements in 3 times:

A:=plottools[polygon]([[0, 1], [2, 1], [2, -1], [0, -1]], style=line,thickness = 2, color=red, numpoints=5000):

B:=plottools[polygon]([[0.03, 1], [2.03, 1.06], [1.97, -0.94], [-0.03, -1]], style=line, thickness = 2, color=blue, numpoints=5000):

plots[display](A,B); 

Manual solution the system using Maple:

Sys:={x+y=1, x+z=2, y+z=3}:

Eq1:=`+`(op(Sys));  # Addition of all three equations

Eq2:=map(t->t/2, Eq1);  # Dividing Eq1 by 2

Eq2-Sys[3];  # Finding  x

Eq2-Sys[2];  # Finding  y

Eq2-Sys[1];  # Finding  z

x:=n->1/2*(x(n-1)+3/x(n-1)):

x(1):=1:

Example - the first 7 terms of the sequense:

seq(x(n), n=1..7);

a:=2046:

a=2*``(a/2);

       2046 = 2*(1023)

If you need to get all the ways of expansion  2046  into 2 factors, then see here

In  Int  in the denominator should be  sqrt(2*g*(a - y(u))) .

I got the final result  1.791629230

What do you mean by an arbitrary partition? If the interval of integration   a..b  divided into  n  equal parts, and the value of the function is taken in the middle of each interval, the Riemann sum is

RiemannSum:=proc(expr, a, b, n)

local x, f, h;

x:=op(indets(expr));

f:=unapply(expr, x);

h:=(b-a)/n;

sum(f(a-h/2+h*i)*h, i=1..n);

end proc:

Examples:

RiemannSum(t^2, 0, 1, 10);

limit(RiemannSum(x^2, 0, 1, n), n=infinity);

For easy use rewrote the program as a procedure. Formal parameters: a  and  b  - sizes of the table,  T - the  duration of the animation,  P - the coordinates of the initial point ,  V - the velocity of the ball, N - (optional) the number of the frames (by default  N=100).

The code of the procedure:

restart;

Billiard:=proc(a::positive, b::positive, T::positive, P::list, V::list, N::posint:=100)

local A, B;

if (P[1]<=0 or P[1]>=a) or (P[2]<=0 or P[2]>=b) then error "Should be 0<=P[1]<=a and 0<=P[2]<=b" fi;

A:=plots[animate](plot,[[2*a/Pi*abs(arcsin(sin(Pi/2/a*(P[1]+V[1]*t)))), 2*b/Pi*abs(arcsin(sin(Pi/2/b*(P[2]+V[2]*t)))), t=0..tau], thickness=5, color=red, numpoints=1000], tau=0..T, frames=N):

B:=plottools[rectangle]([0,b], [a,0], thickness=2, color=green):

plots[display](A, B, scaling=constrained, axes=none)

end proc: 

Examples 1.

Billiard(2, 1, 10, [0.5,0.5], [0.5,-sqrt(3)/2]);

Example 2.  I wonder in what proportions of the parameters the trajectory is closed?

Billiard(4, 2, 40, [1,1], [3/4,2]); 

Here is a simple example of an animation path of a billiard ball on a rectangular table:

restart;

A:=plots[animate](plot,[[3/Pi*abs(arcsin(sin(Pi/2*(0.5+0.2*t)))), 2/Pi*abs(arcsin(sin(Pi/2*(0.5+sqrt(3)/5*t)))), t=0..a], thickness=5, color=red, numpoints=1000], a=0..25, frames=100):

B:=plottools[rectangle]([0,1], [3/2,0], thickness=2, color=green):

plots[display](A, B, scaling=constrained, axes=none);