Rouben Rostamian

MaplePrimes Activity


These are answers submitted by Rouben Rostamian

restart;

P := (u,v) -> -u^4 + 88*u^3*v - 146*u^2*v^2 + 88*u*v^3 - v^4 + 2*u^2 + 40*u*v + 2*v^2 - 1;

proc (u, v) options operator, arrow; -u^4+88*u^3*v-146*u^2*v^2+88*u*v^3-v^4+2*u^2+40*v*u+2*v^2-1 end proc

u0, v0 := I*sqrt(2)/6, -I*sqrt(2)/6;

((1/6)*I)*2^(1/2), -((1/6)*I)*2^(1/2)

Verify that (u0,v0) is on the curve

P(u0,v0);

0

The direct calculation of the slope at (u0, v0) results in division by zero

implicitdiff(P(u,v), v, u);
subs(u=u0,v=v0, %);

(-u^3+66*u^2*v-73*u*v^2+22*v^3+u+10*v)/(-22*u^3+73*u^2*v-66*u*v^2+v^3-10*u-v)

Error, numeric exception: division by zero

So we calculate the slope at (u0, v) and let v approach v0.

implicitdiff(P(u,v), v, u);
subs(u=u0, %):
limit(%, v=v0);

(-u^3+66*u^2*v-73*u*v^2+22*v^3+u+10*v)/(-22*u^3+73*u^2*v-66*u*v^2+v^3-10*u-v)

49/113

The direct evaluation of the second derivative at (u0,v0) also results in division by zero.

So we attempt to evaluate it through a limiting process.

implicitdiff(P(u,v), v, u, u);
subs(u=u0, %):
limit(%, v=v0, left), limit(%, v=v0, right);

(1379*u^8-124432*u^7*v+474796*u^6*v^2-787248*u^5*v^3+898450*u^4*v^4-787248*u^3*v^5+474796*u^2*v^6-124432*u*v^7+1379*v^8-1833*u^6-82844*u^5*v+258133*u^4*v^2-307320*u^3*v^3+258133*u^2*v^4-82844*u*v^5-1833*v^6-1495*u^4-16024*u^3*v+46434*u^2*v^2-16024*u*v^3-1495*v^4-99*u^2-1980*u*v-99*v^2)/(-10648*u^9+105996*u^8*v-447546*u^7*v^2+1026445*u^6*v^3-1352274*u^5*v^4+978663*u^4*v^5-316470*u^3*v^6+13287*u^2*v^7-198*u*v^8+v^9-14520*u^7+94908*u^6*v-237354*u^5*v^2+265701*u^4*v^3-106020*u^3*v^4-9546*u^2*v^5+366*u*v^6-3*v^7-6600*u^5+20580*u^4*v-15486*u^3*v^2-3441*u^2*v^3-138*u*v^4+3*v^5-1000*u^3-300*u^2*v-30*u*v^2-v^3)

infinity, -infinity

The left and right limits of the second derivative at (u0,v0) are unequal.

We conclude that the second derivative is undefined at (u0,v0).

 

Download mw.mw

Comment added later: Some or all of the calculations shown here may be wrong.  See my comment to Kitonum's reply.

Let's say you solve a pde in the unknown u(x,t), as in

pdsol := pdsolve(whatebver, numeric);

To plot the graph of u(0.5,t), we do

pdsol:-plot(x=0.5, t=0..1);

Aside: It will make it easier on others who respond to your question if you include a worksheet along with your question.

I assume that you know the difference between a membrane and a plate.  [A membrane does not support bending moments while a plate does.]

There is not much that Maple can do to help you derive the equations of membranes and plates.  There are various mathematical models that vary depending on the assumptions that you put in.

The simplest plate model is that of Kirchhoff and Love.  See the Wikipedia page for the details of how that model is derived.  Here is the resulting equation in Maple notation, where h is the plate's thickness and w(x,y,t) the the displacement at the point (x,y) at time t, and q(x,y,t) is the applied load:

restart;

The dynamical equation of a Kirchhoff--Love plate

KL__pde := B*(diff(w(x,y,t),x$4)
    + 2*diff(w(x,y,t),x,x,y,y) + diff(w(x,y,t),y$4)) + q(x,y,t) = -2*rho*h*diff(w(x,y,t),t,t);

B*(diff(diff(diff(diff(w(x, y, t), x), x), x), x)+2*(diff(diff(diff(diff(w(x, y, t), x), x), y), y))+diff(diff(diff(diff(w(x, y, t), y), y), y), y))+q(x, y, t) = -2*rho*h*(diff(diff(w(x, y, t), t), t))

where B is

B := 2/3*h^3*E/(1-nu^2);

(2/3)*h^3*E/(-nu^2+1)

In a vibrating plate we have w(x,y,t) = a(x,y)*exp(I*omega*t), where a(x,y) is the amplitude:

subs(q(x,y,t)=0, w(x,y,t)=a(x,y)*exp(I*omega*t), KL__pde):
simplify(%/exp(I*omega*t)/B);

diff(diff(diff(diff(a(x, y), x), x), x), x)+2*(diff(diff(diff(diff(a(x, y), x), x), y), y))+diff(diff(diff(diff(a(x, y), y), y), y), y) = (-3*nu^2+3)*rho*omega^2*a(x, y)/(h^2*E)

Download mw.mw

The eigenvalues ω of that PDE are the vibration frequencies, and the eigenfunctions a(x,y) are the corresponding mode shapes.

As in vv's diagram, replace each of the three radii by a sine curve.

restart;

with(plots):

with(plottools):

f := (x,a) -> a*sin(2*Pi*x);

proc (x, a) options operator, arrow; a*sin(2*Pi*x) end proc

int(sqrt(1+D[1](f)(x,a)^2),x=0..1):
L := simplify(%) assuming a > 0;

2*EllipticE(2*Pi*a/(4*Pi^2*a^2+1)^(1/2))*(4*Pi^2*a^2+1)^(1/2)/Pi

my_a := fsolve(L=2*Pi/3);

.4408865648

C := circle(1, color="Green"):
S1 := plot(f(x,my_a), x=0..1, color="Red"):
S2 := rotate(S1, 2*Pi/3):
S3 := rotate(S1, 4*Pi/3):

display(C, S1, S2, S3, axes=none);

 

Download zz.mw

Maple correctly distinguishes between the n=4 and n≠4 cases because when n=4, the numerator and demominator of the summand in the Fourier series expansion of the solution are zero.

When calculating the Fourier coefficients by hand, one plugs in n=4 before carrying out the required integration.  Maple attempts to "cheat" by calculating the coefficients for general n, and then taking the limit of the general term by letting n→4. That's where things go wrong because n→4 makes no sense when n is an integer.

To see what happens, consider 

(-1 + (-1)^n) / (n-4);
limit(%, n=4);

Maple returns iπ.  The reason for that may be seen more transparently by looking at the stripped-down case

(-1 + (-1)^x)/x;
limit(%, x=0);

where Maple again returns iπ.  That makes sense because (-1)^x as x goes to zero makes no sense without bringing in complex numbers.

Here is the contents of the attached worksheet.  Ignore the artifacts

introduced through this website's interface.  The worksheet itself

is free of those artifacts.

 

The solutions of the first equation are obtained by finding

the intersection points of the graph of m, which is a straight line,

and the graph of tanh(6*m/T) which is a so-called sigmoid.

It is easy to see the graphs intersect at m = 0 regardless of

the value of T.  We conclude that
     m = 0, G = -T*ln(2)

is a solution to your system for all chices of T.

If T >= 6, then m = 0 is the only intersection, and

therefore there is no solution other than the one

shown above when T >= 6.

 

It T < 6, then there are two additional intersections

and we have solutions corresponding

to m+`&+-`(a), for some 0 < a and a < 1 that depends on T.

The proc doit() defined below prints the solution of the system

for any T.   If T < 6,  only one of the three possible

solutions are printed.  

 

restart;

doit := proc(T)
        local eq1, eq2, m, G;
        if not type (T, realcons) then
                return 'procname'(args);
        elif T >= 6 then
                return T, {G=-T*ln(2.0), m=0};
        else
                eq1 := m = tanh(6*m/T);
                eq2 := G = -T*ln(2.0*cosh(6*m/T)) + 3*m^2;
                fsolve({eq1,eq2}, {m=1e-9..1, G=-10..10});
                return T, %;
        end if;
end proc:

for T from 0.1 to 4 by 0.1 do
        doit(T);
end do;

.1, {G = -3., m = 1.}

.2, {G = -3., m = 1.000000000}

.3, {G = -3., m = 1.000000000}

.4, {G = -3., m = 1.000000000}

.5, {G = -3.000000000, m = .9999999999}

.6, {G = -3.000000001, m = .9999999959}

.7, {G = -3.000000027, m = .9999999282}

.8, {G = -3.000000245, m = .9999993882}

.9, {G = -3.000001457, m = .9999967607}

1.0, {G = -3.000006145, m = .9999877098}

1.1, {G = -3.000020125, m = .9999634037}

1.2, {G = -3.000054503, m = .9999091217}

1.3, {G = -3.000127480, m = .9998037082}

1.4, {G = -3.000265626, m = .9996199523}

1.5, {G = -3.000504466, m = .9993256730}

1.6, {G = -3.000888388, m = .9988851653}

1.7, {G = -3.001469979, m = .9982606937}

1.8, {G = -3.002309007, m = .9974138169}

1.9, {G = -3.003471204, m = .9963064344}

2.0, {G = -3.005026998, m = .9949015285}

2.1, {G = -3.007050282, m = .9931636244}

2.2, {G = -3.009617306, m = .9910590175}

2.3, {G = -3.012805698, m = .9885558189}

2.4, {G = -3.016693635, m = .9856238716}

2.5, {G = -3.021359165, m = .9822345760}

2.6, {G = -3.026879655, m = .9783606567}

2.7, {G = -3.033331356, m = .9739758941}

2.8, {G = -3.040789081, m = .9690548347}

2.9, {G = -3.049325946, m = .9635724890}

3.0, {G = -3.059013204, m = .9575040241}

3.1, {G = -3.069920116, m = .9508244480}

3.2, {G = -3.082113883, m = .9435082882}

3.3, {G = -3.095659601, m = .9355292577}

3.4, {G = -3.110620264, m = .9268599051}

3.5, {G = -3.127056751, m = .9174712382}

3.6, {G = -3.145027890, m = .9073323167}

3.7, {G = -3.164590465, m = .8964097942}

3.8, {G = -3.185799289, m = .8846674025}

3.9, {G = -3.208707256, m = .8720653538}

4.0, {G = -3.233365398, m = .8585596366}

for T from 4.5 to 5.0 by 0.05 do
        doit(T);
end do;

4.5, {G = -3.384572858, m = .7755163136}

4.55, {G = -3.402385182, m = .7655763275}

4.60, {G = -3.420707790, m = .7552930332}

4.65, {G = -3.439545617, m = .7446517446}

4.70, {G = -3.458903519, m = .7336364324}

4.75, {G = -3.478786265, m = .7222295433}

4.80, {G = -3.499198538, m = .7104117835}

4.85, {G = -3.520144951, m = .6981618703}

4.90, {G = -3.541630036, m = .6854562243}

4.95, {G = -3.563658247, m = .6722686035}

5.00, {G = -3.586233970, m = .6585696604}

for T from 7.0 to 8.0 by 0.25 do
        doit(T);
end do;

7.0, {G = -4.852030264, m = 0}

7.25, {G = -5.025317059, m = 0}

7.50, {G = -5.198603854, m = 0}

7.75, {G = -5.371890650, m = 0}

8.00, {G = -5.545177445, m = 0}

Download me.mw

The hiker draws on the map an isosceles triangle FAS where FA=SA and the angle at the vertex A is 2*alpha.  Of the two possible such triangles, he picks the one where the vertex A lies closer to him  He draws a circle centered at A and radius AF, and then discards the part of the circle that lies on the other side of the segment FS.  At any point M on the remaining arc, the angle FMS equals alpha.

He repeats the construction on the line segment SK and angle beta.  The intersection P of the two circular arcs constructed this way mark his location since the angles FPS and SPK equal alpha and beta, respectively

Then he reads off the coordinates of the point P on the map in the same way that he read the coordinates of F, S, and K.

if (l1-l2) . (l1-l2) = 0 then
    print("Equal Vectors");
end if;

restart;

eqns :=
        c - d = d - b,
        d = a,
        a^2 + c^2 = 19^2,
        1/2*a*(c-b) = 60;

c-d = d-b, d = a, a^2+c^2 = 361, (1/2)*a*(c-b) = 60

sol := solve({eqns});

{a = RootOf(2*_Z^4-241*_Z^2+3600), b = (1/30)*RootOf(2*_Z^4-241*_Z^2+3600)^3-(181/60)*RootOf(2*_Z^4-241*_Z^2+3600), c = -(1/30)*RootOf(2*_Z^4-241*_Z^2+3600)^3+(301/60)*RootOf(2*_Z^4-241*_Z^2+3600), d = RootOf(2*_Z^4-241*_Z^2+3600)}

subs(sol, sqrt(a^2+b^2)):
simplify(%);

11

 

Download mw.mw

 

Here is a proof.  I suspect that it can be done in a more elegant and more insightful way but I haven't spent much time on it.

restart;

In the diagram below, the red curve is the graph of an arbitrary smooth function f(x).

The green line intersects the curve at A and B.  The region trapped between the

curve and the green line is colored yellow.  

 

Tangent lines to the curve drawn at the points A and B intersect at C, forming

a triangle ABC.  Let
R = (area of the yellow region) / (area of the triangle ABC).

The goal of this worksheet is to show that R converges to 2/3 as B approaches A.

 

Let A=(a,f(a)), B=(b,f(b)).
Equations of tangent lines at x=a and x=b to the curve y=f(x).

tan__1 := y = f(a) + D(f)(a)*(x-a);
tan__2 := y = f(b) + D(f)(b)*(x-b);

y = f(a)+(D(f))(a)*(x-a)

y = f(b)+(D(f))(b)*(x-b)

The intersection point of the tangent lines

inter := simplify(solve({tan__1, tan__2}, {x,y}));

{x = ((D(f))(a)*a-(D(f))(b)*b-f(a)+f(b))/((D(f))(a)-(D(f))(b)), y = (((a-b)*(D(f))(b)+f(b))*(D(f))(a)-f(a)*(D(f))(b))/((D(f))(a)-(D(f))(b))}

Let (X,Y) be the coordinates of the intersection point:

X := simplify(subs(inter, x));
Y := simplify(subs(inter, y));

((D(f))(a)*a-(D(f))(b)*b-f(a)+f(b))/((D(f))(a)-(D(f))(b))

(((a-b)*(D(f))(b)+f(b))*(D(f))(a)-f(a)*(D(f))(b))/((D(f))(a)-(D(f))(b))

The area of the triangle with vertices (a,f(a)), (b,f(b)), (X,Y):

<a,X,b;f(a),Y,f(b);1,1,1>;
A__triangle := simplify(1/2*LinearAlgebra:-Determinant(%));

Matrix(3, 3, {(1, 1) = a, (1, 2) = ((D(f))(a)*a-(D(f))(b)*b-f(a)+f(b))/((D(f))(a)-(D(f))(b)), (1, 3) = b, (2, 1) = f(a), (2, 2) = (((a-b)*(D(f))(b)+f(b))*(D(f))(a)-f(a)*(D(f))(b))/((D(f))(a)-(D(f))(b)), (2, 3) = f(b), (3, 1) = 1, (3, 2) = 1, (3, 3) = 1})

((a-b)*(D(f))(b)-f(a)+f(b))*((a-b)*(D(f))(a)-f(a)+f(b))/(2*(D(f))(a)-2*(D(f))(b))

The equation of the line connecting A and B:

line__ab := y = f(a) + (f(b)-f(a))/(b-a)*(x-a);

y = f(a)+(f(b)-f(a))*(x-a)/(b-a)

The area of the of the region bounded by the line above and the curve y=f(x)
(i.e., the area of the yellow region in the figure)

A__region := int(rhs(line__ab) - f(x), x=a..b);

int(f(a)+(f(b)-f(a))*(x-a)/(b-a)-f(x), x = a .. b)

The ratio of the areas

R := A__region / A__triangle;

(int(f(a)+(f(b)-f(a))*(x-a)/(b-a)-f(x), x = a .. b))*(2*(D(f))(a)-2*(D(f))(b))/(((a-b)*(D(f))(b)-f(a)+f(b))*((a-b)*(D(f))(a)-f(a)+f(b)))

Both the numerator and denominator of R vanish when b=a.  We apply l'Hospital's

rule to find the limit of R as b approaches a:

simplify(diff(A__region, b)):         # derivative of the numerator
simplify(diff(A__triangle, b)):       # derivative of the denominator
simplify(%%/%):                       # l'Hospital's rule
limit(%, b=a);                        # limit of the ratio

2/3

 

 

Download vollers-theorem.mw

 

restart;

with(plots):

p1 := implicitplot([x = y^2 - 3, x = y - y^2], x=-4..1, y=-2..3,
        color=["Red", "Blue"], thickness=3);

p2 := inequal({x > y^2 - 3, x < y - y^2}, x=-4..1, y=-2..3,
        color="DarkSeaGreen", nolines);

display(p1,p2);

 

Download mw.mw

 

I have no idea what "sweep frequency analysis" may mean but I see that your worksheet suffers from bad choices of parameters. Any worksheet that has numbers like 0.0001 and 4000 occurring in the same expression should raise a red flag.  Do you really need those numbers?

I have reproduced your worksheet with a milder set of numbers, and made a few adjustments. See if this is something like what you expect. If so, then change the parameters just a little at a time to make them more to your liking. If not, say what it is that you don't like. A statement like "the two images do not correspond to each other" is not very informative. You need to be more clear about what you mean.

restart;

de := diff(theta(t), t$2) + g*sin(theta(t))/l;

diff(diff(theta(t), t), t)+g*sin(theta(t))/l

g := 1;
T := 40;

1

40

l1 := 0.01 + t/T;
l2 := (1 + 0.01) - t/T;

0.1e-1+(1/40)*t

1.01-(1/40)*t

De1 := subs(l=l1, de);
De2 := subs(l=l2, de);

diff(diff(theta(t), t), t)+sin(theta(t))/(0.1e-1+(1/40)*t)

diff(diff(theta(t), t), t)+sin(theta(t))/(1.01-(1/40)*t)

ICs := theta(0) = 0.1, D(theta)(0) = 0;

theta(0) = .1, (D(theta))(0) = 0

dsol1 := dsolve({De1,ICs}, numeric, output=operator);
dsol2 := dsolve({De2,ICs}, numeric, output=operator);