Rouben Rostamian

A few remarks...

Answered: Rouben Rostamian 8501

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 8501

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 8501

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

Solution...

Answered: Rouben Rostamian 8501

July 07 2020

1 2

>

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 . 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 where the matrix and the right-hand side
(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));

(1)

Let us verify that this reduces to the classic in the

case of a cantilever equation which corresponds to .

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

(2)

Nope, that didn't work!

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

order and therefore plugging into it is not the right things to do. Instead,

we divide the equation by , and then take the limit as 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 .

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

in Figure 2 which corresponds to the choice of . 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

>

Download mw.mw

Finding the angle α between P0P1 a...

Answered: Rouben Rostamian 8501

July 06 2020

0 1

Finding the angle α between P₀P₁ and P₀P₂ is a simple matter of applying the law of cosines in the triangle P₀P₁P₂, which is expressed as
||P₂ - P₁||² = ||P₁ - P₀||² + ||P₂ - P₀||² - 2 ||P₁ - P₀|| ||P₂ - P₀|| cos(α),
whence
cos(α) = [ L₁₂² - L₁₀² - L₂₀² ] / [ 2 L₁₀ L₂₀ ].

Is this what you need?...

Answered: Rouben Rostamian 8501

July 05 2020

0 2

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[]]

Maybe this?...

Answered: Rouben Rostamian 8501

July 03 2020

2 3

restart;
n:= 10: alpha:= -60: beta:= 20:
Qs:= alpha + beta*P:
P_val := 4:
printf("There are %a identical firms.  Supply is given by Qs = %a.  When P = %a then quantity supplied is equal to %a.  Calculate the ...\n", n, Qs, P_val, eval(Qs, P=P_val));

Upon execution, this prints:

There are 10 identical firms. Supply is given by Qs = 20*P-60. When P = 4 then quantity supplied is equal to 20. Calculate the ...

Needs better care...

Answered: Rouben Rostamian 8501

July 02 2020

1 0

@666 jvbasha You can begin by paying closer attention to your code. In AN.mw you have

t(x) := sum(p^i*t[i](x), i = 0 .. L);

Do you see that it defines t in terms of t ? Same with the definition of g. I have not read the code any further.

plot(x^2, x=0..1, labels=[xi, `f '`(...

Answered: Rouben Rostamian 8501

June 30 2020

2 1

plot(x^2, x=0..1, labels=[ xi, `f '`(xi) ] );

Corrections...

Answered: Rouben Rostamian 8501

June 29 2020

3 1

Your differential equations are of the first order. That means that you may specify initial conditions of the values of the unknowns, but not their derivatives. The solution will determine what those derivatives are. This is what math says, not maple. Your specification of initial values of derivatives is illegal. Remove them. Keep only ics := {h1(0.0001) = 1, h2(0.0001) = 1, h3(0.0001) = 1};
Your system of equations consisting of eq2, eq3, eq4 is a system of ordinary differential equations (ODEs). Maple's ODE solver is called dsolve. You are applying pdsolve which is designed for solving parial differential equations (PDEs). So change the pdsolve line to
dsol := dsolve(sys1 union ics, numeric);
Then you may plot the solution h1(z) through
plots:-odeplot(dsol, [z, h1(z)], z=0..1);
When you plot h2 and h3, you will see that they are identical to h1. That's because your equations and initial conditions are completely symmetrical with respect to h1, h2, h3, and therefore the h1, h2, and h3 are identical.
Additional note: At the top of your worksheet you have:
with(PDEtools): with(LinearAlgebra): with(DifferentialGeometry): with(VectorCalculus): with(DEtools):
These are entirely irrelevant to the rest of the worksheet. Remove!

Needs a refined initmesh...

Answered: Rouben Rostamian 8501

June 28 2020

2 0

>

restart;

>	kernelopts(version);

>	eq1 := diff(f_n(zeta), zeta$2)-h_n(zeta)(diff(f_n(zeta), zeta))-f_n(zeta)^2+g_n(zeta)^2 = 0; eq2 := diff(g_n(zeta), zeta$2)-h_n(zeta)(diff(g_n(zeta), zeta))-2f_n(zeta)g_n(zeta) = 0; eq3 := 2*f_n(zeta)+diff(h_n(zeta), zeta) = 0; ics := f_n(0) = 0, f_n(L) = 0, g_n(0) = 1, g_n(L) = 0, h_n(0) = 0;

diff(diff(f_n(zeta), zeta), zeta)-h_n(zeta)*(diff(f_n(zeta), zeta))-f_n(zeta)^2+g_n(zeta)^2 = 0

diff(diff(g_n(zeta), zeta), zeta)-h_n(zeta)*(diff(g_n(zeta), zeta))-2*f_n(zeta)*g_n(zeta) = 0

>	L := 10.0; dsol := dsolve({eq1,eq2,eq3,ics}, numeric, initmesh = 16): # default initmesh is 8 plots:-odeplot(dsol, [[zeta,f_n(zeta)], [zeta, g_n(zeta)], [zeta,h_n(zeta)]], color=[red,green,blue]);

Download mw.mw

Answer...

Answered: Rouben Rostamian 8501

June 27 2020

3 22

>

restart;

>	pde := diff(w(x,t), x$4) + diff(w(x,t), t$2)/c^2 = 0;

diff(diff(diff(diff(w(x, t), x), x), x), x)+(diff(diff(w(x, t), t), t))/c^2 = 0

We look for solutions of the form :

>	eval(pde, w(x,t)=u(x)cos(omegat)); de_tmp := expand(%/cos(omega*t));

(diff(diff(diff(diff(u(x), x), x), x), x))*cos(omega*t)-u(x)*omega^2*cos(omega*t)/c^2 = 0

diff(diff(diff(diff(u(x), x), x), x), x)-u(x)*omega^2/c^2 = 0

The coefficient is not convenient. It's better to work with dimensionless quantities.

Therefore we introduce a new unknown, , related to through

>	mu_def := omega = c/L^2*mu^2;

omega = c*mu^2/L^2

Then express the differential equation de_tmp in terms of

>	de := eval(de_tmp, mu_def);

diff(diff(diff(diff(u(x), x), x), x), x)-u(x)*mu^4/L^4 = 0

Here is the general solution of the de:

>	dsolve(de): U := unapply(rhs(%), x);

proc (x) options operator, arrow; _C1*exp(mu*x/L)+_C2*exp(-mu*x/L)+_C3*sin(mu*x/L)+_C4*cos(mu*x/L) end proc

We want to apply the boundary conditions , but of course
Maple will pick the trivial solution . The get around that, we temporarily suspend the
boundary condition and instead substitute a boundary condition , where
is an unspecified symbol. Applying the four boundary conditions, we obtain four equations in

the four unknowns which we solve:

>	U(0) = 0, D(U)(0)=0, (D@@2)(U)(0) = a, U(L) = 0: C_vals := solve({%}, {_C1, _C2, _C3, _C4});

${_C1 = -(1/2)*a*L^2*(-cos(mu)+sin(mu)+exp(-mu))/(mu^2*(exp(mu)-exp(-mu)-2*sin(mu))), _C2 = (1/2)*a*L^2*(exp(mu)-cos(mu)-sin(mu))/((exp(mu)-exp(-mu)-2*sin(mu))*mu^2), _C3 = (1/2)*a*L^2*(exp(mu)-2*cos(mu)+exp(-mu))/(mu^2*(exp(mu)-exp(-mu)-2*sin(mu))), _C4 = -(1/2)*a*L^2/mu^2}$

Plug these into the expression for , compute the derivative, and set the derivative to zero
in order to assert the boundary condition which we had left out.

>	eval(U(x), C_vals): diff(%, x): eval(%, x=L) = 0: simplify(%);

a*L*(exp(2*mu)*cos(mu)-2*exp(mu)+cos(mu))/(mu*(exp(2*mu)-2*sin(mu)*exp(mu)-1)) = 0

The denominator of that equation is a multiple of and therefore nonzero.
We conclude that

>	exp(2mu)cos(mu) - 2*exp(mu) + cos(mu) = 0; eq_tmp := isolate(%, cos(mu));

cos(mu) = 2*exp(mu)/(exp(2*mu)+1)

Let's verify that the right-hand side of this equation is

>	simplify(1/cosh(mu) - rhs(eq_tmp));

Okay, therefore eq_tmp is equivalent to:

>	eq := cos(mu) = 1/cosh(mu);

This equation has infinitely many roots, as can be seen in this graph:

>	plot([cos(mu),1/cosh(mu)], mu=0..5*Pi, color=[blue,red], thickness=3, size=[700,300]);

Here are the first few roots:

>	seq(fsolve(cos(mu)=1/cosh(mu), mu=(2n+1)Pi/2), n=1..5);

Call the roots .. For each root there is a mode function given by

>	subs(C_vals, mu=mu[j], U(x));

-(1/2)*a*L^2*(-cos(mu[j])+sin(mu[j])+exp(-mu[j]))*exp(mu[j]*x/L)/(mu[j]^2*(exp(mu[j])-exp(-mu[j])-2*sin(mu[j])))+(1/2)*a*L^2*(exp(mu[j])-cos(mu[j])-sin(mu[j]))*exp(-mu[j]*x/L)/((exp(mu[j])-exp(-mu[j])-2*sin(mu[j]))*mu[j]^2)+(1/2)*a*L^2*(exp(mu[j])-2*cos(mu[j])+exp(-mu[j]))*sin(mu[j]*x/L)/(mu[j]^2*(exp(mu[j])-exp(-mu[j])-2*sin(mu[j])))-(1/2)*a*L^2*cos(mu[j]*x/L)/mu[j]^2

where the factors may be dropped since a mode function is determined modulo a

multiplicative constant. Also for each there corresponds a frequency given through

the equation mu_def, that is, "`omega__j`=(c*`mu__j`^(2))/(L^(2))." You should be all set now for eigenfunction expansion

of the solutions of your PDE.

>

Download mw.mw

Solution...

Answered: Rouben Rostamian 8501

June 27 2020

4 1

>

restart;

>	d := 6/7; # the delay

Plot the history function

>	p[0] := plot(cos(t), t=-d..0):

The solution on the interval picks up the history function (colored green)

>	diff(x(t),t,t) + 2/25diff(x(t),t) + 4x(t) = (1/25diff(x(t),t,t) + 2/3diff(x(t),t) + 5/2x(t))cos(t-d), x(0)=1, D(x)(0)=0: dsol := dsolve({%}, numeric, output=operator, range=0..d): myx[1] := eval(x, dsol): mydx[1] := eval(D(x), dsol): p[1] := plot(myx[1](t), t=0..d):

On each subsequent interval , the differential equation picks up the solution

calculated in the preceding interval (colored green)

>

for i from 2 to 14 do
sys := diff(x(t),t,t) + 2/25*diff(x(t),t) + 4*x(t)
= (1/25*diff(x(t),t,t) + 2/3*diff(x(t),t) + 5/2*x(t))*myx[i-1](t-d),
x((i-1)*d)=myx[i-1]((i-1)*d), D(x)((i-1)*d)=mydx[i-1]((i-1)*d):
dsol := dsolve({sys}, numeric, output=operator, range=(i-1)*d..i*d):
myx[i] := eval(x, dsol):
mydx[i] := eval(D(x), dsol):
p[i] := plot(myx[i](t), t=(i-1)*d..i*d);
end do:
plots:-display(seq(p[j], j=0..i-1));

Download mw.mw

Solution...

Answered: Rouben Rostamian 8501

June 25 2020

2 4

Here is the solution with no explanations. See if you can figure it out.

>

restart;

>	C1 := x^2 + y^2 = L^2; C2 := (x-R)^2 + y^2 = R^2;

>	solve(C1, y): c1 := %[1];

>	solve(C2, y): c2 := %[1];

>	X := solve(c1=c2, x);

>	A1 := int(c1, x=X..L) assuming L>0, 2*R>L;

-(1/8)*L^2*(4*arcsin((1/2)*L/R)*R^2-2*Pi*R^2+L*(-L^2+4*R^2)^(1/2))/R^2

>	A2 := int(c2, x=0..X) assuming L>0, 2*R>L;

(1/8)*(2*R^4*Pi+4*R^4*arcsin((1/2)*(L^2-2*R^2)/R^2)+L^3*(-L^2+4*R^2)^(1/2)-2*L*(-L^2+4*R^2)^(1/2)*R^2)/R^2

>	A := simplify(2*(A1+A2));

(1/2)*(-2*arcsin((1/2)*L/R)+Pi)*L^2-(1/2)*L*(-L^2+4*R^2)^(1/2)+(arcsin((1/2)*(L^2-2*R^2)/R^2)+(1/2)*Pi)*R^2

>	A = Pi*R^2/2: eval(%, R=1): fsolve(%, L=0..2);

Does not exist...

Answered: Rouben Rostamian 8501

June 23 2020

0 1

Such a rational function does not exist! Is that a trick question?

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

8501 Reputation

23 Badges

MaplePrimes Activity

These are answers submitted by Rouben Rostamian

A few remarks...

Plotting...

Complete solution...

Solution...

Finding the angle α between P0P1 a...

Is this what you need?...

Maybe this?...

Needs better care...

plot(x^2, x=0..1, labels=[xi, `f '`(...

Corrections...

Needs a refined initmesh...

Answer...

Solution...

Solution...

Does not exist...

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