Rouben Rostamian

Animation...

Answered: Rouben Rostamian 7711

September 13 2020

0 0

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.

>

restart;

>	with(geometry):

Define the ellipse through its semi-major axis length and focal length

>	a, c := 3, 2;

The eccentricity

>

e := c/a;

The foci:

>	point(B,0,0): point(C,2*c,0):

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:

Animate:

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

Download mw.mw

Reduction to an ODE...

Answered: Rouben Rostamian 7711

August 26 2020

1 3

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.

>

restart;

>	Do := f -> diff(f, r $ 2) + diff(f, phi $ 2)/r^2 - cos(phi)diff(f, phi)/(r^2sin(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);

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

>	exp(-_c[1]t)sin(phi)r(_C2BesselJ(1,sqrt(_c[1]/nu)sin(phi)r) + _C3BesselY(1,sqrt(_c[1]/nu)sin(phi)r)): psi(t,r,phi) = %w(rsin(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))

Let's get rid of the irrelevant multiplicative factor:

>	expr := r^2nusin(phi)^2 / exp(-_c[1]t) expr_tmp;

Simplify and then introduce :

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

That's a fourth order ODE in the unknown . 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
, etc.

Download mw.mw

Try this...

Answered: Rouben Rostamian 7711

August 14 2020

3 1

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.

>

restart;

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

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

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

>	k := x->piecewise(x<=0.5,3-1x,x>0.5,2.5); c := x->piecewise(x<=0.5,1.2-0.4x,x>0.5,1); rho := x->piecewise(x<=0.5,20-10x,x>0.5,15); f0 := x->t0x; 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

>

pde_sys;

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

Download mw.mw

Contact sales...

Answered: Rouben Rostamian 7711

August 10 2020

1 0

Go to https://www.maplesoft.com/pricing/

Solution...

Answered: Rouben Rostamian 7711

August 03 2020

2 5

>

restart;

I will use instead of for the sake of clarity;

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

>	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 &DifferentialD;y &DifferentialD;x" :

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

Moment about the axis: "`M__y `=(∫)[x=0]^(1) (∫)[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(%);

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

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

The moment of inertia about the origin: "`I__o`=(∫)[x=0]^(1) (∫)[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;

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

>	evalf(I__c);

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

No symbolic solution...

Answered: Rouben Rostamian 7711

July 31 2020

0 0

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:

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.

The general case...

Answered: Rouben Rostamian 7711

July 27 2020

1 2

Here it is...

Answered: Rouben Rostamian 7711

July 25 2020

1 0

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

No calculations necessary...

Answered: Rouben Rostamian 7711

July 25 2020

3 6

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.

Answer...

Answered: Rouben Rostamian 7711

July 22 2020

1 0

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

>	doit(3,4,5);

>	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: 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

Solution...

Answered: Rouben Rostamian 7711

July 18 2020

3 3

An illustration...

Answered: Rouben Rostamian 7711

July 15 2020

2 1

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 in the first equation by the corresponding value in the second
equation, that is 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
, 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

A few remarks...

Answered: Rouben Rostamian 7711

July 15 2020

1 8

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.

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.
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.
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.
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.

Plotting...

Answered: Rouben Rostamian 7711

July 13 2020

3 4

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

Complete solution...

Answered: Rouben Rostamian 7711

July 08 2020

1 7

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 := rhoAdiff(u(x,t),t,t) + EI*diff(u(x,t),x$4) = 0;

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

>	eval(pde, u(x,t)=X(x)cos(omegat)); 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 . In the subsequent calculations we will find . Then we

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

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

functions instead. We let be the beam's deflection in each of the

three spans.

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

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 Here we construct the coefficient of the system's matrix. The
system's right-hand side, which is a vector of zeros, is captured in the vector 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;

_rtable[18446883936173081894]

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

That yields a transcendental equation for the unknown .

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

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;

Plug the span lengths into the characteristic equation

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

>	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

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

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

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

of the beam, that is, the coefficients of the original PDE are immaterial.

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

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

therefore the determinant of 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 ), 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

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post

7711 Reputation

22 Badges

MaplePrimes Activity

These are answers submitted by Rouben Rostamian

Animation...

Reduction to an ODE...

Try this...

Contact sales...

Solution...

No symbolic solution...

The general case...

Here it is...

No calculations necessary...

Answer...

Solution...

An illustration...

A few remarks...

Plotting...

Complete solution...

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Save this setting as your default sorting preference?
Note: You can change your preference any time in your account settings
Don't show this again