Rouben Rostamian

MaplePrimes Activity

These are answers submitted by Rouben Rostamian

Here is a slightly modified version of your worksheet.  The determinant grows too fast for efficient numerical calculations.  Therefore I have scaled it by dividing by cosh(...) which does not affect the roots since cosh is never zero.  Then you see that you can calculate the first five roots without any fuss.



If you look closely at your equations, you will see that at some places you have x, and at some other places you have x(t).  They should all be x(t).  Same goes for y(t) and u(t).

I have fixed these in the attached worksheet.  Note the technique—to define a differential equation dx/dt = F(x,y,u), we begin with defining the function (not an expression!) F, and then express the equation as dx(t)/dt = F(x(t), y(t), u(t)).

But note that the solution to your system is not defined in the entire length of the time interval that you have requested.  Since you have x, y, and u in the denominators, the solution stops when any of these variables reaches zero.



p := 2*u*x*y^2;
q := (x-1)/x;



doth := simplify(-(sqrt(6)*x^3*(y-1)^2+u*((3+sqrt(6))*x^3*(y-1)^2+3*y^2-9*x*y^2+3*x^2*(y^2+2*y-1)))/p);


ddotphi := simplify((y-1)*(sqrt(6)*x*y^2*(x-1)+u*(3*y^2-(6+sqrt(6))*x*y^2+x^2*((3+sqrt(6))*y^2+6*y-3)))/p);


x*(1-x)*doth - x^2*ddotphi/sqrt(6):
F1 := unapply(%, [x,y,u]);    # the right-hand side of ode1

proc (x, y, u) options operator, arrow; (1/2)*(1-x)*(-x^3*(y-1)^2*(u+1)*6^(1/2)-3*(x-1)*((y-1)^2*x^2+2*x*y^2-y^2)*u)/(u*y^2)-(1/12)*x*(y-1)*(x*y^2*(x-1)*(u+1)*6^(1/2)+3*((x-1)^2*y^2+2*x^2*y-x^2)*u)*6^(1/2)/(u*y^2) end proc

ode1 := diff(x(t),t) = F1(x(t),y(t),u(t));

diff(x(t), t) = (1/2)*(1-x(t))*(-x(t)^3*(y(t)-1)^2*(u(t)+1)*6^(1/2)-3*(x(t)-1)*((y(t)-1)^2*x(t)^2+2*x(t)*y(t)^2-y(t)^2)*u(t))/(u(t)*y(t)^2)-(1/12)*x(t)*(y(t)-1)*(x(t)*y(t)^2*(x(t)-1)*(u(t)+1)*6^(1/2)+3*((x(t)-1)^2*y(t)^2+2*x(t)^2*y(t)-x(t)^2)*u(t))*6^(1/2)/(u(t)*y(t)^2)

F2 := unapply(%, [x,y,u]);        # the right-hand side of ode2

proc (x, y, u) options operator, arrow; (1/2)*(y-1)*(-x^3*(y-1)^2*(u+1)*6^(1/2)-3*(x-1)*((y-1)^2*x^2+2*x*y^2-y^2)*u)/(y*u*x) end proc

ode2 := diff(y(t), t) = F2(x(t),y(t),u(t));

diff(y(t), t) = (1/2)*(y(t)-1)*(-x(t)^3*(y(t)-1)^2*(u(t)+1)*6^(1/2)-3*(x(t)-1)*((y(t)-1)^2*x(t)^2+2*x(t)*y(t)^2-y(t)^2)*u(t))/(y(t)*u(t)*x(t))

F3 := unapply(%, [x,y,u]);          # the right-hand side of ode3

proc (x, y, u) options operator, arrow; (1/2)*6^(1/2)*u*(u+1)*(x-1)/x end proc

ode3 := diff(u(t), t) = F3(x(t),y(t),u(t));

diff(u(t), t) = (1/2)*6^(1/2)*u(t)*(u(t)+1)*(x(t)-1)/x(t)


{t, diff(u(t), t), diff(x(t), t), diff(y(t), t), u(t), x(t), y(t)}

DEplot3d([ode1, ode2, ode3], [x(t), y(t), u(t)], t = -2 .. 2, x = 0 .. 1, y = 0 .. 1, u = -1 .. 0, [[x(0) = .9, y(0) = .25, u(0) = -.95], [x(0) = .9, y(0) = .25, u(0) = -.5], [x(0) = .85, y(0) = .2, u(0) = -.15], [x(0) = .25, y(0) = .5, u(0) = -.5], [x(0) = .7, y(0) = .5, u(0) = -.6]], linecolor = [red, green, black, navy, maroon], thickness = 2, axes = boxed, labels = [x(t), y(t), u(t)], stepsize = 0.1e-1, orientation = [30, 65]);

Warning, plot may be incomplete, the following errors(s) were issued:
   cannot evaluate the solution further right of .16790195, probably a singularity
Warning, plot may be incomplete, the following errors(s) were issued:
   cannot evaluate the solution further right of .51702724e-1, probably a singularity
Warning, plot may be incomplete, the following errors(s) were issued:
   cannot evaluate the solution further right of .57085115e-2, probably a singularity
Warning, plot may be incomplete, the following errors(s) were issued:
   cannot evaluate the solution further right of .54099962, probably a singularity




The English term for "exinscrit circle" in the second sentence is  "excircle".

Here is an animation of the excircle tangent to the side AC.  We see that it goes through the ellipse's vertex, as asserted.  I haven't attempted a proof but let me add an extra observation: It appears that the locus of the excircle's center is a vertical line.



Define the ellipse through its semi-major axis length a and focal length "c." 

a, c := 3, 2;

3, 2

The eccentricity

e := c/a;


The foci:


Proc to draw one animation frame

frame := proc(t)
        local r, A, T, E, obj, circ;
        r := a*(1-e^2)/(1 - e*cos(t));  # equation of ellipse in polar coords
        point(A, r*cos(t), r*sin(t));
        triangle(T, [B,C,A]);
        ellipse(E, [foci=[B,C], distance=2*a]);
        excircle(obj, T);
        circ := %[1];
        draw([T,E,circ], color=[red,"Green",blue], printtext=true, axes=none);
end proc:


nframes := 30:
plots:-display([seq(frame(i/nframes*Pi/2.0), i=1..nframes)], insequence, scaling=constrained);








I don't like the idea of using both psi and psi[2] in the same Maple worksheet.  In this particular case that's harmless but often it is not.  So I have changed psi[2] to w.


Do := f -> diff(f, r $ 2) + diff(f, phi $ 2)/r^2 - cos(phi)*diff(f, phi)/(r^2*sin(phi));

proc (f) options operator, arrow; diff(f, `$`(r, 2))+(diff(f, `$`(phi, 2)))/r^2-cos(phi)*(diff(f, phi))/(r^2*sin(phi)) end proc

pde := Do(Do(psi(t, r, phi)) - diff(psi(t, r, phi), t)/nu);

diff(diff(diff(diff(psi(t, r, phi), r), r), r), r)+(diff(diff(diff(diff(psi(t, r, phi), phi), phi), r), r))/r^2-4*(diff(diff(diff(psi(t, r, phi), phi), phi), r))/r^3+6*(diff(diff(psi(t, r, phi), phi), phi))/r^4-cos(phi)*(diff(diff(diff(psi(t, r, phi), phi), r), r))/(r^2*sin(phi))+4*cos(phi)*(diff(diff(psi(t, r, phi), phi), r))/(r^3*sin(phi))-6*cos(phi)*(diff(psi(t, r, phi), phi))/(r^4*sin(phi))-(diff(diff(diff(psi(t, r, phi), r), r), t))/nu+(diff(diff(diff(diff(psi(t, r, phi), phi), phi), r), r)+(diff(diff(diff(diff(psi(t, r, phi), phi), phi), phi), phi))/r^2+2*(diff(diff(psi(t, r, phi), phi), phi))/r^2-cos(phi)*(diff(diff(diff(psi(t, r, phi), phi), phi), phi))/(r^2*sin(phi))+2*cos(phi)^2*(diff(diff(psi(t, r, phi), phi), phi))/(r^2*sin(phi)^2)-2*cos(phi)*(diff(psi(t, r, phi), phi))/(r^2*sin(phi))-2*cos(phi)^3*(diff(psi(t, r, phi), phi))/(r^2*sin(phi)^3)-(diff(diff(diff(psi(t, r, phi), phi), phi), t))/nu)/r^2-cos(phi)*(diff(diff(diff(psi(t, r, phi), phi), r), r)+(diff(diff(diff(psi(t, r, phi), phi), phi), phi))/r^2+(diff(psi(t, r, phi), phi))/r^2-cos(phi)*(diff(diff(psi(t, r, phi), phi), phi))/(r^2*sin(phi))+cos(phi)^2*(diff(psi(t, r, phi), phi))/(r^2*sin(phi)^2)-(diff(diff(psi(t, r, phi), phi), t))/nu)/(r^2*sin(phi))

Substitute a special form of psi into the PDE.  This introduces a new unknown w:

exp(-_c[1]*t)*sin(phi)*r*(_C2*BesselJ(1,sqrt(_c[1]/nu)*sin(phi)*r) + _C3*BesselY(1,sqrt(_c[1]/nu)*sin(phi)*r)):
psi(t,r,phi) = %*w(r*sin(phi));
expr_tmp := simplify(eval(pde, %));

psi(t, r, phi) = exp(-_c[1]*t)*sin(phi)*r*(_C2*BesselJ(1, (_c[1]/nu)^(1/2)*sin(phi)*r)+_C3*BesselY(1, (_c[1]/nu)^(1/2)*sin(phi)*r))*w(r*sin(phi))

2*exp(-_c[1]*t)*((5/2)*_C2*((-(1/5)*cos(phi)^2*r^2*_c[1]+(1/5)*r^2*_c[1]-(3/5)*nu)*(D(w))(r*sin(phi))+((2/5)*nu*r*(cos(phi)-1)*(cos(phi)+1)*((D@@3)(w))(r*sin(phi))+((cos(phi)^2*r^2*_c[1]-r^2*_c[1]+(3/5)*nu)*((D@@2)(w))(r*sin(phi))-(1/5)*((D@@4)(w))(r*sin(phi))*nu*r^2*(cos(phi)-1)*(cos(phi)+1))*sin(phi))*r)*BesselJ(1, (_c[1]/nu)^(1/2)*sin(phi)*r)+(5/2)*_C3*((-(1/5)*cos(phi)^2*r^2*_c[1]+(1/5)*r^2*_c[1]-(3/5)*nu)*(D(w))(r*sin(phi))+((2/5)*nu*r*(cos(phi)-1)*(cos(phi)+1)*((D@@3)(w))(r*sin(phi))+((cos(phi)^2*r^2*_c[1]-r^2*_c[1]+(3/5)*nu)*((D@@2)(w))(r*sin(phi))-(1/5)*((D@@4)(w))(r*sin(phi))*nu*r^2*(cos(phi)-1)*(cos(phi)+1))*sin(phi))*r)*BesselY(1, (_c[1]/nu)^(1/2)*sin(phi)*r)+r^3*(_c[1]/nu)^(1/2)*sin(phi)*(cos(phi)-1)*(cos(phi)+1)*(BesselJ(0, (_c[1]/nu)^(1/2)*sin(phi)*r)*_C2+BesselY(0, (_c[1]/nu)^(1/2)*sin(phi)*r)*_C3)*((D(w))(r*sin(phi))*_c[1]-2*((D@@3)(w))(r*sin(phi))*nu))/(r^2*nu*sin(phi)^2)

Let's get rid of the irrelevant multiplicative factor:

expr := r^2*nu*sin(phi)^2 / exp(-_c[1]*t) * expr_tmp;

5*_C2*((-(1/5)*cos(phi)^2*r^2*_c[1]+(1/5)*r^2*_c[1]-(3/5)*nu)*(D(w))(r*sin(phi))+((2/5)*nu*r*(cos(phi)-1)*(cos(phi)+1)*((D@@3)(w))(r*sin(phi))+((cos(phi)^2*r^2*_c[1]-r^2*_c[1]+(3/5)*nu)*((D@@2)(w))(r*sin(phi))-(1/5)*((D@@4)(w))(r*sin(phi))*nu*r^2*(cos(phi)-1)*(cos(phi)+1))*sin(phi))*r)*BesselJ(1, (_c[1]/nu)^(1/2)*sin(phi)*r)+5*_C3*((-(1/5)*cos(phi)^2*r^2*_c[1]+(1/5)*r^2*_c[1]-(3/5)*nu)*(D(w))(r*sin(phi))+((2/5)*nu*r*(cos(phi)-1)*(cos(phi)+1)*((D@@3)(w))(r*sin(phi))+((cos(phi)^2*r^2*_c[1]-r^2*_c[1]+(3/5)*nu)*((D@@2)(w))(r*sin(phi))-(1/5)*((D@@4)(w))(r*sin(phi))*nu*r^2*(cos(phi)-1)*(cos(phi)+1))*sin(phi))*r)*BesselY(1, (_c[1]/nu)^(1/2)*sin(phi)*r)+2*r^3*(_c[1]/nu)^(1/2)*sin(phi)*(cos(phi)-1)*(cos(phi)+1)*(BesselJ(0, (_c[1]/nu)^(1/2)*sin(phi)*r)*_C2+BesselY(0, (_c[1]/nu)^(1/2)*sin(phi)*r)*_C3)*((D(w))(r*sin(phi))*_c[1]-2*((D@@3)(w))(r*sin(phi))*nu)

Simplify and then introduce eta = r*sin(phi):

simplify(expr, {cos(phi)^2 = 1 - sin(phi)^2}):
algsubs(r*sin(phi)=eta, %):
ode := collect(%, [(D@@4)(w)(eta), (D@@3)(w)(eta), (D@@2)(w)(eta), (D@@2)(w)(eta)], simplify);

eta^3*nu*(BesselJ(1, (_c[1]/nu)^(1/2)*eta)*_C2+BesselY(1, (_c[1]/nu)^(1/2)*eta)*_C3)*((D@@4)(w))(eta)-2*(BesselJ(1, (_c[1]/nu)^(1/2)*eta)*_C2+BesselY(1, (_c[1]/nu)^(1/2)*eta)*_C3-2*eta*(_c[1]/nu)^(1/2)*(BesselJ(0, (_c[1]/nu)^(1/2)*eta)*_C2+BesselY(0, (_c[1]/nu)^(1/2)*eta)*_C3))*eta^2*nu*((D@@3)(w))(eta)+3*eta*(-(5/3)*_c[1]*eta^2+nu)*(BesselJ(1, (_c[1]/nu)^(1/2)*eta)*_C2+BesselY(1, (_c[1]/nu)^(1/2)*eta)*_C3)*((D@@2)(w))(eta)-3*(_C2*(-(1/3)*_c[1]*eta^2+nu)*BesselJ(1, (_c[1]/nu)^(1/2)*eta)+_C3*(-(1/3)*_c[1]*eta^2+nu)*BesselY(1, (_c[1]/nu)^(1/2)*eta)+(2/3)*(_c[1]/nu)^(1/2)*eta^3*_c[1]*(BesselJ(0, (_c[1]/nu)^(1/2)*eta)*_C2+BesselY(0, (_c[1]/nu)^(1/2)*eta)*_C3))*(D(w))(eta)

That's a fourth order ODE in the unknown w(eta).  Maple is unable to obtain the general solution of

that ODE.  Perhaps solutions can be obtained if you specify the values of some of the constants
_c__1, _C__3, etc.




Look at this post and the ensuing discussion on ways to handle discontinuous coefficients in the heat equation.   Converting the heat equation into a system of two first order equations (suggested in that discussion by vv) works in your case.  Here it is.


pde_sys :=
        k(x)*diff(T(x,t),x) = Q(x,t),
        diff(Q(x,t),x) = c(x)*rho(x)*diff(T(x,t),t);

k(x)*(diff(T(x, t), x)) = Q(x, t), diff(Q(x, t), x) = c(x)*rho(x)*(diff(T(x, t), t))

bc := D[1](T)(0,t) = q0, T(1,t) = t0;

(D[1](T))(0, t) = q0, T(1, t) = t0

ic := T(x,0) = f0(x);

T(x, 0) = f0(x)

k        := x->piecewise(x<=0.5,3-1*x,x>0.5,2.5);
c        := x->piecewise(x<=0.5,1.2-0.4*x,x>0.5,1);
rho      := x->piecewise(x<=0.5,20-10*x,x>0.5,15);
f0       := x->t0*x;
t0       := 20;
q0       := -2;

k := proc (x) options operator, arrow; piecewise(x <= .5, 3-x, .5 < x, 2.5) end proc

c := proc (x) options operator, arrow; piecewise(x <= .5, 1.2+(-1)*.4*x, .5 < x, 1) end proc

proc (x) options operator, arrow; piecewise(x <= .5, 20-10*x, .5 < x, 15) end proc

proc (x) options operator, arrow; t0*x end proc




piecewise(x <= .5, 3-x, .5 < x, 2.5)*(diff(T(x, t), x)) = Q(x, t), diff(Q(x, t), x) = piecewise(x <= .5, 1.2-.4*x, .5 < x, 1)*piecewise(x <= .5, 20-10*x, .5 < x, 15)*(diff(T(x, t), t))

pdsol := pdsolve({pde_sys}, {ic, bc}, numeric, spacestep=0.01);


tmax := 20;


pdsol:-animate(T(x,t), t=0..tmax, frames=100);

pdsol:-animate(Q(x,t), t=0..tmax, frames=100);








I will use x instead of t for the sake of clarity;

f := x -> sqrt(x^(2/3) - x^2);

proc (x) options operator, arrow; sqrt(x^(2/3)-x^2) end proc

plot([f(x),-f(x)], x=0..1, scaling=constrained);

The area is given by the double integral  "(&int;)[x=0  ]^(1)(&int;)[y=-f(x)]^(f(x))1 &DifferentialD;y &DifferentialD;x":

area := int(int(1, y=-f(x)..f(x)), x=0..1);


Moment about the y axis: "`M__y `=(&int;)[x=0]^(1)  (&int;)[y=-f(x)]^(f(x))x &DifferentialD;y &DifferentialD;x"

int(int(x, y=-f(x)..f(x)), x=0..1):
convert(%, hypergeom):
M__y := simplify(%);


The x coordinate of the center of mass:

x__c := M__y / area;




The moment of inertia about the origin:  "`I__o`=(&int;)[x=0]^(1)  (&int;)[y=-f(x)]^(f(x))((x^(2)+y^(2) )&DifferentialD;y &DifferentialD;x"

I__o := int(int((x^2+y^2), y=-f(x)..f(x)), x=0..1);


The moment of intertia about the center of mass obtained by applying the parallel axis theorem:

I__c := I__o - area * x__c^2;






Note added later:

Going from the Cartesian (x,y) coordinates to the polar (r,t) coordinates where x=r*cos(t), y=r*sin(t),  the equation of the lamina's boundary changes over to r = sqrt(cos(t)), where t goes from −π/2 to π/2. All the previous calculations may be repeated just as easily with double integrals in polar coordinates.



From the way you have phrased your question, it appears that you are under the impression that Maple is not printing the solution because it is too long.  But that's not true.  Maple is not printing a solution because your equation has no solution expressible in symbolic form.  The following much shorter example shows a similar outcome:


de := diff(y(x),x$4) + (a*x + b)*y(x) = 1;

diff(diff(diff(diff(y(x), x), x), x), x)+(a*x+b)*y(x) = 1

bc := y(0)=0, D(y)(0)=0, y(1)=0, D(y)(1)=0;

y(0) = 0, (D(y))(0) = 0, y(1) = 0, (D(y))(1) = 0

ans := dsolve({de,bc});

In such situations often the only recourse is to look for a numerical approximation to the solution.  For that, you need to begin with specifying the numerical values of all the constants that enter your equation.


Let's say "{`phi__n` : n=1,2...}" is an orthonormal basis in some Hilbert space H.

Then any `in`(f, H) may be written as f = sum((f, `&phi;__n`)*`&phi;__n`, n = 1 .. infinity).  That's as general as it gets.

If you have a particular set of basis functions, just replace `&phi;__n` with those.

If you replace them with sine or cosine you will get the traditional Fourier series.


If you have something else in mind, then you will need to be more specific.


A := < 4, 3, -5 >:
B := < 3, -1, 4 >:
plots:-arrow(A, B-A, scaling=constrained, view=[-5..5, -5..5, -5..5]);

The ratio of the areas of the rectangles EFGH and SJIK equals (ML/NL)^2.
This remains true even if you replace the rectangle EFGH by any convex shape.

The proc doit( ) does what you have asked.  It calculates the value of the cell (i,j) in the nxn case.


doit := proc(i,j,n)
        local upper_left, d;
        upper_left := proc(i,j,n)
                local d, z;
                d := i + j - 1;
                z := ( d^2+1+(d-1)*(-1)^d )/2;
                `if`(type(d, odd), z + (j-1),  z - (j-1));
        end proc:
        d := i + j - 1;
        if d <= n then
                % + 2*(d-n)*n - (d-n)^2;
        end if;
end proc:

Examples of use:





doit(30, 40, 100);


Formatted listing of the 5x5 case:

n := 5:
for i from 1 to n do
        for j from 1 to n do
                printf("%4d", doit(i, j, n)):
        end do:
end do:
unassign('n', 'i', 'j');

   1   2   6   7  15
   3   5   8  14  16
   4   9  13  17  22
  10  12  18  21  23
  11  19  20  24  25

Formatted listing of the 8x8 case:

n := 8:
for i from 1 to n do
        for j from 1 to n do
                printf("%4d", doit(i, j, n)):
        end do:
end do:
unassign('n', 'i', 'j');

   1   2   6   7  15  16  28  29
   3   5   8  14  17  27  30  43
   4   9  13  18  26  31  42  44
  10  12  19  25  32  41  45  54
  11  20  24  33  40  46  53  55
  21  23  34  39  47  52  56  61
  22  35  38  48  51  57  60  62
  36  37  49  50  58  59  63  64





ans := FourierSeries(f(x), x = -Pi .. Pi, p, series = sine);

Sum((int(f(x)*sin((1/2)*i*(x+Pi)), x = -Pi .. Pi))*sin((1/2)*i*(x+Pi))/Pi, i = 1 .. p)

subs(i=n, ans);

Sum((int(f(x)*sin((1/2)*n*(x+Pi)), x = -Pi .. Pi))*sin((1/2)*n*(x+Pi))/Pi, n = 1 .. p)


The attached worksheed illustrates how to handle one of your several cases.  The other cases should work in the same way.



A vertical cylider as a parametric surface:

cyl1 := plot3d([cos(t), sin(t), z], t=0..2*Pi, z=-2..2, color=red):

A horizontal cylider as a parametric surface:

cyl2 := plot3d([cos(t), y, sin(t)], t=0..2*Pi, y=-2..2, color=blue):

The intersection curve lies  on both cylinders, therefore when we are on the curve,

we may replace z in the first equation by the corresponding value in the second
equation, that is "sin(t)."  We get

curve1 := tubeplot([cos(t), sin(t), sin(t)], t=0..2*Pi, radius=0.06, color=yellow):

But a surface's parametrization is not unique.  We could have parametrized cyl1 as
[cos(t), sin(t), -z], in which curse the intersection curve would have been:

curve2 := tubeplot([cos(t), sin(t), -sin(t)], t=0..2*Pi, radius=0.06, color=yellow):

Here is all four objects put together:

display(cyl1, cyl2, curve1, curve2, scaling=constrained, style=surface);




As Carl has pointed out, the message regarding output limit is harmless.  That's provided that it were correct.  The problem is that it is not correct.  See if you can fix it following these hints.

  1. A symbolic computational program such as Maple is quite different from purely numerical programming languages such as C, Fortran, Matlab.  In those programming languages you enter the values of constants, such as E, I, kt, kr, etc, right at the beginning of the program.  In a symbolic computational program it is always better to delay for as long as possible the entering of the numerical values.
    In your worksheet you have entered numerical values for kt and kr at the very beginning.  That's bad.  Just use the symbols kt and kr for them until their numerical values become absolutely necessary.  That will fix the output limit exceeded issue.
  2. The interface conditions at the supports are missing the E.I  factors.  Need to insert those and define values for E and I.  But as I have noted above, do not insert the values at the beginning.  Delay for as long as possible.
  3. Note that the symbol I in maple stands for the square root of -1.  We don't want to redefine that symbol.  Instead, note that the parameters E and I always enter the equations as the product E*I.  Therefore it is best to introduce a single two-letter symbol EI that stands for that product.  Carry out the symbolic calculations with EI, and substitute a numerical value for it at the end.
  4. If you are following the formulas in the appendix of the paper by Akeson, et. al., for the interface conditions, beware that those formulas are incorrect.  For instance, subtracting the third and fifth equations on page 89 yields d(w1)/dz = 0, which is not true.  (A beam's slope is almost never zero at a support.)  Since this is the subject of your doctoral research, I suggest that you derive those equations from basic principles for yourself rather than replying on that paper.

There may be some other issues with the worksheet, but see if you can address these ones first.


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