Rouben Rostamian

MaplePrimes Activity


These are answers submitted by Rouben Rostamian

restart;

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  "(∫)[x=0  ]^(1)(∫)[y=-f(x)]^(f(x))1 ⅆy ⅆx":

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

1

Moment about the y axis: "`M__y `=(∫)[x=0]^(1)  (∫)[y=-f(x)]^(f(x))x ⅆy ⅆx"

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

(2/5)*GAMMA(3/4)^2*2^(1/2)/Pi^(1/2)

The x coordinate of the center of mass:

x__c := M__y / area;

(2/5)*GAMMA(3/4)^2*2^(1/2)/Pi^(1/2)

evalf(x__c);

.4792560936

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

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

(1/8)*Pi

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

I__c := I__o - area * x__c^2;

-(8/25)*GAMMA(3/4)^4/Pi+(1/8)*Pi

evalf(I__c);

.1630126787

 

 

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.

Download mw.mw

 

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:

restart;

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, `φ__n`)*`φ__n`, n = 1 .. infinity).  That's as general as it gets.

If you have a particular set of basis functions, just replace `φ__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.

restart;

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
                upper_left(i,j);
        else
                upper_left(n+1-j,n+1-i,n);
                % + 2*(d-n)*n - (d-n)^2;
        end if;
end proc:

Examples of use:

doit(3,2,5);

9

doit(3,4,5);

17

doit(30, 40, 100);

2386

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:
        printf("\n");
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:
        printf("\n");
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
 

Download mw.mw

 

restart;

with(OrthogonalExpansions):

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.

restart;

with(plots):

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);

 

Download mw.mw

 

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.

 

In your calculation you arrive at
eq := (some expression in mu, L[1], L[2], L[3]) = 0;

After setting L[1]=0, you get an equation in mu, L[2], L[3].

To plot, introduce the variables

rho = L[3]/L[2]              # the length ratios
lambda = mu*L[2]         # scaled eigenvalue

Then the equation reduces to an expression only in mu and rho
which then we may plot with implictplot.

eval(eq, L[1]=0);
eval(%, L[3]=rho*L[2]);
eval(%, mu=lambda/L[2]);
plots:-implicitplot(%, rho=0..5, lambda=0..5);

Here is the calculation that reproduces the graphs in Figure 3.22(a) on page 70.

Note that the boundary conditions specified on page 88 seem to be incorrect (or perhaps I am misinterpreting them.)  For that reason I have applied the boundary conditions in the way that make sense to me.  I suggest that you use these ones instead unless you figure out the correct interpretation of those on page 88.

Modal analysis of a 3-span beam

Here we formulate and solve for the modal shapes of transverse vibration of a 3-span

free-pinned-pinned-free beam.  The graphs of the first three modes agree with those

in Figure 3.22 on page 70 of the 2007 article of
Henrik Åkesson, Tatiana Smirnova, Thomas Lagö, and Lars Håkansson.

 

Note: Their boundary conditions, given in (5.1) on page 88, don't look correct to me,

so here I do it my way.

 

Rouben Rostamian

2020-07-08

restart;

#Digits := 15:

with(LinearAlgebra):

pde := rho*A*diff(u(x,t),t,t) + EI*diff(u(x,t),x$4) = 0;

rho*A*(diff(diff(u(x, t), t), t))+EI*(diff(diff(diff(diff(u(x, t), x), x), x), x)) = 0

To find the fundamental modes, look for a time-harmonic solution:

eval(pde, u(x,t)=X(x)*cos(omega*t));
ode_tmp := isolate(%, diff(X(x),x$4));

-rho*A*X(x)*omega^2*cos(omega*t)+EI*(diff(diff(diff(diff(X(x), x), x), x), x))*cos(omega*t) = 0

diff(diff(diff(diff(X(x), x), x), x), x) = rho*A*X(x)*omega^2/EI

So our task has been reduced to solving the ODE above.  The coefficient on the right-hand side is

unnecessarily messy, so we name it mu^4.  In the subsequent calculations we will find mu.  Then we

can find "omega," if needed, from the relationship below:

coeff(rhs(ode_tmp), X(x)) = mu^4;
omega_vs_mu := isolate(%, omega^2);

rho*A*omega^2/EI = mu^4

omega^2 = mu^4*EI/(rho*A)

ode := eval(ode_tmp, omega_vs_mu);

diff(diff(diff(diff(X(x), x), x), x), x) = X(x)*mu^4

Solve the ODE:

dsolve(ode);

X(x) = _C1*exp(-mu*x)+_C2*exp(mu*x)+_C3*sin(mu*x)+_C4*cos(mu*x)

Maple presents the solution in terms of exponential functions.  We prefer hyperbolic

functions instead. We let X__1(x), X__2(x), X__3(x) be the beam's deflection in each of the

three spans.

X[1] := x -> a[1]*sin(mu*x) + b[1]*cos(mu*x) + c[1]*sinh(mu*x) + d[1]*cosh(mu*x);
X[2] := x -> a[2]*sin(mu*x) + b[2]*cos(mu*x) + c[2]*sinh(mu*x) + d[2]*cosh(mu*x);
X[3] := x -> a[3]*sin(mu*x) + b[3]*cos(mu*x) + c[3]*sinh(mu*x) + d[3]*cosh(mu*x);

proc (x) options operator, arrow; a[1]*sin(mu*x)+b[1]*cos(mu*x)+c[1]*sinh(mu*x)+d[1]*cosh(mu*x) end proc

proc (x) options operator, arrow; a[2]*sin(mu*x)+b[2]*cos(mu*x)+c[2]*sinh(mu*x)+d[2]*cosh(mu*x) end proc

proc (x) options operator, arrow; a[3]*sin(mu*x)+b[3]*cos(mu*x)+c[3]*sinh(mu*x)+d[3]*cosh(mu*x) end proc

so here is the overall deflection:

modal_shape := piecewise(x<L[1], X[1](x), x<L[1]+L[2], X[2](x), X[3](x));

modal_shape := piecewise(x < L[1], a[1]*sin(mu*x)+b[1]*cos(mu*x)+c[1]*sinh(mu*x)+d[1]*cosh(mu*x), x < L[1]+L[2], a[2]*sin(mu*x)+b[2]*cos(mu*x)+c[2]*sinh(mu*x)+d[2]*cosh(mu*x), a[3]*sin(mu*x)+b[3]*cos(mu*x)+c[3]*sinh(mu*x)+d[3]*cosh(mu*x))

We have three spans, there are two boundary conditions at the end of each span,

and therefore we have 12 boundary conditions altogether.  These are almost,

but not quite, the same as those given on page 88 of the paper cited above.

bcs :=
(D@@2)(X[1])(0) = 0,   # zero moment at x=0
(D@@3)(X[1])(0) = 0,   # zero shear  at x=0

X[1](L[1]) = 0,        # zero displacement at x=L[1]
X[2](L[1]) = 0,        # zero displacement at x=L[1]
D(X[1])(L[1]) = D(X[2])(L[1]),            # continuous slope at x=L[1]
(D@@2)(X[1])(L[1]) = (D@@2)(X[2])(L[1]),  # no external moment at x=L[1]

X[2](L[1]+L[2]) = 0,        # zero displacement at x=L[1]+L[2]
X[3](L[1]+L[2]) = 0,        # zero displacement at x=L[1]+L[2]
D(X[2])(L[1]+L[2]) = D(X[3])(L[1]+L[2]),            # continuous slope at x=L[1]+L[2]
(D@@2)(X[2])(L[1]+L[2]) = (D@@2)(X[3])(L[1]+L[2]),  # no external moment at x=L[1]+L[2]

(D@@2)(X[3])(L[1]+L[2]+L[3]) = 0,   # zero moment at x=L[1]+L[2]+L[3]
(D@@3)(X[3])(L[1]+L[2]+L[3]) = 0:   # zero shear  at x=L[1]+L[2]+L[3]

The 12 conditions given above constitute a system of 12 linear, homogeneous equations in the

12 variables "`a__1`,`b__1`. ..., `c__3`, `d__3`."  Here we construct the coefficient M of the system's matrix.  The
system's right-hand side, which is a vector of zeros, is captured in the vector MM which is not

used in the subsequent calculations.

M, MM := LinearAlgebra:-GenerateMatrix({bcs},
  [ a[1], b[1], c[1], d[1], a[2], b[2], c[2], d[2], a[3], b[3], c[3], d[3] ]):

Here is the 12x12 matrix M:

M;

_rtable[18446883936173081894]

We set the determinant of M to zero in order to get a nontrivial solution to our homogeneous system.

That yields a transcendental equation for the unknown mu.

Determinant(M):
expand(%):
simplify(%):
characteristic_equation := %/(8*mu^16);  # get rid of the unnecessary factor of 8*mu^16

(((-2*sinh(mu*L[2])*sin(mu*L[2])*cos(mu*L[3])+sin(mu*L[3])*(cosh(mu*L[2])*sin(mu*L[2])-sinh(mu*L[2])*cos(mu*L[2])))*cos(mu*L[1])+sin(mu*L[1])*((cosh(mu*L[2])*sin(mu*L[2])-sinh(mu*L[2])*cos(mu*L[2]))*cos(mu*L[3])+sin(mu*L[3])*(cosh(mu*L[2])*cos(mu*L[2])-1)))*cosh(mu*L[3])+(-sinh(mu*L[3])*(cosh(mu*L[2])*sin(mu*L[2])-sinh(mu*L[2])*cos(mu*L[2]))*cos(mu*L[3])-2*sinh(mu*L[2])*sin(mu*L[2]))*cos(mu*L[1])-sin(mu*L[1])*(sinh(mu*L[3])*(cosh(mu*L[2])*cos(mu*L[2])-1)*cos(mu*L[3])-cosh(mu*L[2])*sin(mu*L[2])+sinh(mu*L[2])*cos(mu*L[2])))*cosh(mu*L[1])+(-((cosh(mu*L[2])*sin(mu*L[2])-sinh(mu*L[2])*cos(mu*L[2]))*cos(mu*L[3])+sin(mu*L[3])*(cosh(mu*L[2])*cos(mu*L[2])-1))*sinh(mu*L[1])*cos(mu*L[1])-2*sinh(mu*L[2])*sin(mu*L[2])*cos(mu*L[3])+sin(mu*L[3])*(cosh(mu*L[2])*sin(mu*L[2])-sinh(mu*L[2])*cos(mu*L[2])))*cosh(mu*L[3])+(sinh(mu*L[3])*(cosh(mu*L[2])*cos(mu*L[2])-1)*cos(mu*L[3])-cosh(mu*L[2])*sin(mu*L[2])+sinh(mu*L[2])*cos(mu*L[2]))*sinh(mu*L[1])*cos(mu*L[1])-sinh(mu*L[3])*(cosh(mu*L[2])*sin(mu*L[2])-sinh(mu*L[2])*cos(mu*L[2]))*cos(mu*L[3])-2*sinh(mu*L[2])*sin(mu*L[2])

We are ready now for some number crunching.  For the span lengths we take the values

given in the caption of Figure 2.12 on page 28.

span_lengths := L[1]=3.5, L[2]=5.0, L[3]=21.5;

L[1] = 3.5, L[2] = 5.0, L[3] = 21.5

Plug the span lengths into the characteristic equation

characteristic_equation_numeric := eval(characteristic_equation, {span_lengths});

(((-2*sinh(5.0*mu)*sin(5.0*mu)*cos(21.5*mu)+sin(21.5*mu)*(cosh(5.0*mu)*sin(5.0*mu)-sinh(5.0*mu)*cos(5.0*mu)))*cos(3.5*mu)+sin(3.5*mu)*((cosh(5.0*mu)*sin(5.0*mu)-sinh(5.0*mu)*cos(5.0*mu))*cos(21.5*mu)+sin(21.5*mu)*(cosh(5.0*mu)*cos(5.0*mu)-1)))*cosh(21.5*mu)+(-sinh(21.5*mu)*(cosh(5.0*mu)*sin(5.0*mu)-sinh(5.0*mu)*cos(5.0*mu))*cos(21.5*mu)-2*sinh(5.0*mu)*sin(5.0*mu))*cos(3.5*mu)-sin(3.5*mu)*(sinh(21.5*mu)*(cosh(5.0*mu)*cos(5.0*mu)-1)*cos(21.5*mu)-cosh(5.0*mu)*sin(5.0*mu)+sinh(5.0*mu)*cos(5.0*mu)))*cosh(3.5*mu)+(-((cosh(5.0*mu)*sin(5.0*mu)-sinh(5.0*mu)*cos(5.0*mu))*cos(21.5*mu)+sin(21.5*mu)*(cosh(5.0*mu)*cos(5.0*mu)-1))*sinh(3.5*mu)*cos(3.5*mu)-2*sinh(5.0*mu)*sin(5.0*mu)*cos(21.5*mu)+sin(21.5*mu)*(cosh(5.0*mu)*sin(5.0*mu)-sinh(5.0*mu)*cos(5.0*mu)))*cosh(21.5*mu)+(sinh(21.5*mu)*(cosh(5.0*mu)*cos(5.0*mu)-1)*cos(21.5*mu)-cosh(5.0*mu)*sin(5.0*mu)+sinh(5.0*mu)*cos(5.0*mu))*sinh(3.5*mu)*cos(3.5*mu)-sinh(21.5*mu)*(cosh(5.0*mu)*sin(5.0*mu)-sinh(5.0*mu)*cos(5.0*mu))*cos(21.5*mu)-2*sinh(5.0*mu)*sin(5.0*mu)

plot(characteristic_equation_numeric, mu=0..0.5);

It is difficult to see the roots in the graph above.  Therefore we limit the vertical range:

plot(characteristic_equation_numeric, mu=0..0.5, view=-4..4);

We proceed to find the first few roots.

 

Important!  The roots depends on the choice of the span lengths "`L__1`, `L__2`, `L__3`."

If these are changed, then the initial guesses in the calculations below should be
re-adjusted according to what is seen in the graph above:

mu__1 := fsolve(characteristic_equation_numeric, mu=0.08);
mu__2 := fsolve(characteristic_equation_numeric, mu=0.20);
mu__3 := fsolve(characteristic_equation_numeric, mu=0.35);
mu__4 := fsolve(characteristic_equation_numeric, mu=0.42);

0.8148236435e-1

.2065743153

.3465175842

.4130255340

We write a proc that calculates and plots the modal shape for any of the mu values

calculated above.  Note that the shapes don't depend on the physical properties

of the beam, that is, the coefficients rho, E, I, A of the original PDE are immaterial.

 

If the mu values shown above were calculated perfectly, then the determinant of M

would be exactly zero.  In reality, the values of mu are only approximations, and

therefore the determinant of M is not exactly zero and depending on the values

of the span lengths, it may be quite far from zero.  The proc prints out the value

of the determinant.  If it is not small (like less that 10^(-5)),  then increase the value

of Digits and try again.

plot_mode := proc(mu)
  local coeff_vals, total_length, my_M, my_mode_shape;
  total_length := add(rhs~({span_lengths}));
  my_M := eval(M, {span_lengths, :-mu=mu});
  printf("Determinant is %g.  Increase Digists if not near zero!\n", Determinant(my_M));
  NullSpace(my_M)[1]:
  <a[1], b[1], c[1], d[1], a[2], b[2], c[2], d[2], a[3], b[3], c[3], d[3]> =~ %:
  coeff_vals := convert(%, set)[];
  #print(coeff_vals);
  my_mode_shape := eval(modal_shape, {:-mu=mu, coeff_vals, span_lengths}):
  eval(my_mode_shape, x=total_length);
  # the (-%) below normalizes the graph so that its value is -1 at the right end
  plot(my_mode_shape/(-%), x=0..total_length);
end proc:

p1 := plot_mode(mu__1);

Determinant is -1.88829e-26.  Increase Digists if not near zero!

p2 := plot_mode(mu__2);

Determinant is -6.39035e-15.  Increase Digists if not near zero!

p3 := plot_mode(mu__3);

Determinant is -2.99856e-07.  Increase Digists if not near zero!

Here are the first three modal plotted together.  This agrees with Figure 3.22(a) in the cited article.

plots:-display([p1,p2,p3], color=["Red","Green","Blue"],
   legend=[mode1, mode2, mode3], axes=boxed);

 

 

 

Download beam-modal-shapes.mw

 

 

restart;

with(LinearAlgebra):

It is clear from the calculations presented in the paper that the quantities

that matter are mu = beta*L, nu = a/L, rho = m__attached/m__beam.  Their equations would have

been much neater if written in terms of mu, nu, rho.  That's what I will do in what
follows and I suggest that you do it this way in your dissertation.   

 

Here is the system of linear equations (14) through (21):

sys :=
a[2]+a[4]=0,
a[1]+a[3]=0,
-b[1]*sin(mu)-b[2]*cos(mu)+b[3]*sinh(mu)+b[4]*cosh(mu)=0,
b[1]*(-cos(mu) + mu*rho*sin(mu)) + b[2]*(sin(mu)+mu*rho*cos(mu))
  + b[3]*(cosh(mu)+mu*rho*sinh(mu)) + b[4]*(sinh(mu)+mu*rho*cosh(mu)) = 0,
a[1]*sin(mu*nu) + a[2]*cos(mu*nu) + a[3]*sinh(mu*nu) + a[4]*cosh(mu*nu) = 0,
b[1]*sin(mu*nu) + b[2]*cos(mu*nu) + b[3]*sinh(mu*nu) + b[4]*cosh(mu*nu) = 0,
a[1]*cos(mu*nu) - a[2]*sin(mu*nu) + a[3]*cosh(mu*nu) + a[4]*sinh(mu*nu)
  - b[1]*cos(mu*nu) + b[2]*sin(mu*nu) - b[3]*cosh(mu*nu) - b[4]*sinh(mu*nu) = 0,
-a[1]*sin(mu*nu) - a[2]*cos(mu*nu) + a[3]*sinh(mu*nu) + a[4]*cosh(mu*nu)
  + b[1]*sin(mu*nu) + b[2]*cos(mu*nu) - b[3]*sinh(mu*nu) - b[4]*cosh(mu*nu) = 0:

We write this in matrix form "M*X=B," where the matrix M and the right-hand side B 
(which is a vector of zeros) are computed through:

M,B :=GenerateMatrix({sys}, [a[1],a[2],a[3],a[4],b[1],b[2],b[3],b[4]]):

To have nonzero solutions to this homogeneous system, we need the determinant

to be zero. That gives us the equation that determines the system's eigenvalues:

eq := simplify(Determinant(M));

(-4*sinh(mu)*sin(mu)*mu*rho-4*cosh(mu)*(sin(mu)+mu*rho*cos(mu)))*cosh(mu*nu)^2+(((8*mu*rho*sin(mu)-4*cos(mu))*sinh(mu)+4*cosh(mu)*sin(mu))*cos(mu*nu)+((4*mu*rho*cos(mu)+4*sin(mu))*sinh(mu)+4*mu*rho*sin(mu)*cosh(mu))*sinh(mu*nu)-8*((mu*rho*cos(mu)+(1/2)*sin(mu))*sinh(mu)+(1/2)*cosh(mu)*cos(mu)+1/2)*sin(mu*nu))*cosh(mu*nu)+((-4*mu*rho*sin(mu)+4*cos(mu))*sinh(mu)+4*mu*rho*cos(mu)*cosh(mu))*cos(mu*nu)^2+((-8*mu*rho*sin(mu)*cosh(mu)+4*cosh(mu)*cos(mu)-4*sinh(mu)*sin(mu)+4)*sinh(mu*nu)+4*((sin(mu)+mu*rho*cos(mu))*sinh(mu)+mu*rho*sin(mu)*cosh(mu))*sin(mu*nu))*cos(mu*nu)+8*sinh(mu*nu)*((1/2)*sinh(mu)*cos(mu)+cosh(mu)*(mu*rho*cos(mu)+(1/2)*sin(mu)))*sin(mu*nu)

(1)

Let us verify that this reduces to the classic cos(mu)*cosh(mu)+1 = 0 in the

case of a cantilever equation which corresponds to rho = 0, nu = 0.

eval(eq, {rho=0, nu=0});

0

(2)

Nope, that didn't work!

 

Looking for the source of trouble, we see that the equation comes with a factor of

order nu^3and therefore plugging nu = 0 into it is not the right things to do.  Instead,

we divide the equation by nu^3, and then take the limit as nu goes to zero:

eval(eq, rho=0):
limit(%/nu^3, nu=0):
factor(%);

-(8/3)*mu^3*(cosh(mu)*cos(mu)+1)

(3)

and there we have the expected equation cos(mu)*cosh(mu)+1 = 0.

 

Now with that puzzle out of the way, we turn to producing the plot

in Figure 2 which corresponds to the choice of rho = .2.  Here it is.
You don't need a for-loop for that.

eval(eq, rho=0.2):
plots:-implicitplot(%, nu=0..1, mu=0..18);

 

You can obtain Figures 3 through 6 in the same way by varying the value of "rho."

 

 

Download mw.mw

 

Finding the angle α between P0P1 and P0P2 is a simple matter of applying the law of cosines in the triangle P0P1P2, which is expressed as
||P2 - P1||2 = ||P1 - P0||2 + ||P2 - P0||2 - 2 ||P1 - P0|| ||P2 - P0|| cos(α),
whence
cos(α) = [ L122 - L102 - L202 ] / [ 2 L10 L20 ].

You have f both on the right-hand side and on the left-hand side of the definition of f.  I suppose that on the right-hand side you meant to have e instead of f.  I have fixed that.  The list L is the desired result.

restart;
f:=a*u[]*u[1]*u[2]+b*u[1]*u[2]+c*v[]*v[1]*u[2]+d*u[]*v[1]*v[2]+d*u[]*v[1]*u[2]+e*u[]*v[];
    f := a u[] u[1] u[2] + c u[2] v[] v[1] + d u[] u[2] v[1] + d u[] v[1] v[2] + b u[1] u[2] + e u[] v[]
collect(f,{u[1],u[2],v[1],v[2]},distributed);
op([1..-1], %);
L := eval([%], [u[1]=1,u[2]=1,v[1]=1,v[2]=1]);
       (a u[] + b) u[1] u[2] + (c v[] + d u[]) u[2] v[1] + d u[] v[1] v[2] + e u[] v[]
       (a u[] + b) u[1] u[2], (c v[] + d u[]) u[2] v[1],  d u[] v[1] v[2], e u[] v[]
       L := [a u[] + b, c v[] + d u[], d u[], e u[] v[]]

 

5 6 7 8 9 10 11 Last Page 7 of 44