Rouben Rostamian

MaplePrimes Activity


These are answers submitted by Rouben Rostamian

vv has already provided a solution. Here is an alternative approach.

restart;

f:=x->1/(x^2-2*x);

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

p1 := plot(f(x), x=-3..5, discont, view=-3..3, scaling=constrained, color=blue);

f is symmetric about the line x = 1

is( f(x+1) = f(-x+1) );

true

Therefore, the center of the circle lies on the line x = 1. The circle touches

the graph at 1, -1, and consequently, the y coordinate of its center is R-1,
where R is the radius. We conclude that the circle's equation is

C := (x-1)^2 + (y-R+1)^2 = R^2;

(x-1)^2+(y-R+1)^2 = R^2

We calculate the points of intersection of the circle and the graph of f:

xy := solve({C, y=f(x)}, [x,y]);

[[x = 1, y = -1], [x = RootOf(1+_Z^4-4*_Z^3+(-2*R+5)*_Z^2+(4*R-2)*_Z), y = -RootOf(1+_Z^4-4*_Z^3+(-2*R+5)*_Z^2+(4*R-2)*_Z)^2+2*R+2*RootOf(1+_Z^4-4*_Z^3+(-2*R+5)*_Z^2+(4*R-2)*_Z)-1]]

Let us isolate the argument of the RootOf:

select(has, xy, RootOf):
eval(x, %[]):
P := op(1, %);

1+_Z^4-4*_Z^3+(-2*R+5)*_Z^2+(4*R-2)*_Z

The quartic polynomial P has 4 roots in general, but when the

circle is tangent to the graph of f,  those 4 roots collapse

into 2 roots.  That happens when the discriminant of the quartic

is zero.

discrim(P, _Z) = 0;
solve(%);

512*R^5-768*R^4-768*R^3+640*R^2+672*R+144 = 0

3/2, 3/2, -1/2, -1/2, -1/2

We conclude that the radius is 3/2 and therefore the y coordinate of its center is 1/2.

Here is its graph

plots:-display(p1, plottools:-circle([1,1/2], 3/2, color=red));


 

 

 

Download mw.mw

 

In your worksheet you have
motioneq  = (3.2)

That "=" should be ":=".

I know nothing about MapleSim, so I cannot help you in that respect.

 

 

 

Unless you have a very ancient version of maple, this should do what you want:

int~(<sin(x), cos(x)>, x=0..2);

I can think of a dozen or two different answers to your extremely vague question and it's likely that they will be unsatisfactory as long as you refuse to show your code. Anyway, here is one possibility:

restart;

E := phi -> 1 / (1+phi^2);

proc (phi) options operator, arrow; 1/(1+phi^2) end proc

E(0.4);

.8620689655

plot(E(phi), phi=0..1);

 

 

 

 

Download mw.mw

Based on the information that you have provided, this is all that we can do:

restart;

de := diff(x(t),t) = A(t)*x(t) + B(t)*u(t);

diff(x(t), t) = A(t)*x(t)+B(t)*u(t)

dsol := dsolve(de, x(t));

x(t) = (Int(B(t)*u(t)*exp(-(Int(A(t), t))), t)+_C1)*exp(Int(A(t), t))

y(t) = eval( C(t)*x(t) + E(t)*u(t), dsol);

y(t) = C(t)*(Int(B(t)*u(t)*exp(-(Int(A(t), t))), t)+_C1)*exp(Int(A(t), t))+E(t)*u(t)

If that's not what you want, then you will have to be more specific about what you are looking for.

 

restart;

x := [x__1,y__1,theta__1,x__2,y__2,theta__2];

[x__1, y__1, theta__1, x__2, y__2, theta__2]

f := [cos(theta__1),sin(theta__1),0,0,0,0];

[cos(theta__1), sin(theta__1), 0, 0, 0, 0]

J := Matrix(6, 6, (i,j) -> diff(f[i], x[j]));

Matrix(6, 6, {(1, 1) = 0, (1, 2) = 0, (1, 3) = -sin(`#msub(mi("&theta;",fontstyle = "normal"),mi("1"))`), (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = cos(`#msub(mi("&theta;",fontstyle = "normal"),mi("1"))`), (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0, (4, 5) = 0, (4, 6) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = 0, (5, 5) = 0, (5, 6) = 0, (6, 1) = 0, (6, 2) = 0, (6, 3) = 0, (6, 4) = 0, (6, 5) = 0, (6, 6) = 0})

 

Download mw.mw

 

Here are two different ways to achieve that.

Method 1: Use the EulerLagrange function in the VariationalCalculus package:

VariationalCalculus:-EulerLagrange(L, t, [ x(t), y(t), alpha(t), phi1(t), phi2(t) ] );
remove(type, %, equation);

Method 2: Load the Physics package:

with(Physics):

That redefines diff so that your code should work unless you have other errors in it.

There are other ways as well but the two that I have shown here are the quickest. 

In the future, when you ask a question, include your actual Maple worksheet.  Showing just a picture, as you have done, is not very helpful since I cannot execute that code to test my suggestions before replying.

 

The polygonplot3d help page describes three different ways of invoking that function.  Their descriptions are somewhat difficult to decipher, so here is an expanded description.

In the first method, we pass the vertex coordinates to polygonplot3d as a single list of lists, each internal list providing the coordinates of one vertex of the polygon:

plots:-polygonplot3d( [ [x1,y1,z1], [x2,y2,z2], ..., [xn,yn,zn] ]);

In the second method we pass an n by 3 matrix to polygonplot3d, were the rows of the matrix are the coordinates of the polygon's vertices.

In the third method we pass three vectors of length n each to polygonplot3d.  If these vectors are called U, V, W, then ( U[k], V[k], W[k] ) is taken to be the coordinates of the vertex number k.

plots:-polygonplot3d(U,V,W);

Your call to polygonplot3d invokes the thrid method. Thus, the first plotted point has coordinates [17/2,-17/2,0] as you have observed.  That's not what you want.

What you want is to invoke the first method which requires the vertices to be given as lists, not vectors.  There are two ways you may approach this: either provide vertices as lists, as in

A1 := [17/2, 0, 0];
A2 := [-17/2, 0, 0];
B1 := [0, 11/2, sqrt(166/4)];
plots:-polygonplot3d( [A1,A2,B1]);

or provide the vertices as vectors as you already have, and have maple convert them to lists before plotting, as in

A1 := Vector([17/2, 0, 0]);
A2 := Vector([-17/2, 0, 0]);
B1 := Vector([0, 11/2, sqrt(166/4)]);
plots:-polygonplot3d(convert~([A1,A2,B1], list));

 

 

For the weak formulation, you multiply the equation by a so-called "test function", let's say v, and then integrate the various terms by parts in order to reduce the differentiation order of the u function and increase the differentiation order of the v function.  For instance the fourth order derivative term, after integration by parts, changes to ∇2 u . ∇2v .  I assume that you know how to do this, otherwise you will need to study some very basic ideas related to weak solutions of PDEs.  The formal computation above makes sense if the functions u and v are in the Sobolev space H2

After you have done that, you introduce the Galerkin approximation of functions in Hwhereby your PDE reduces to a system of linear ODEs in the time variable.

For the computational implementation you will need to pick basis functions for the Galerkin approximation.  That depends on the dimension of the space variable.  Is your space one or two or three dimensional? I am guessing that you are interested in the two-dimensional case.  In that case we will want triangular elements with Argyris shape functions.  (Do a search for "Argyris element" in google.)  Argyris shape functions are designed to connect across element boundaries as C1 functions. That ensures that the Galerkin approximation forms a proper subspace of H2.

 

  • What is the unknown?  Is it u or w?
  • What is f? Don't you mean f(x,y)?
  • This being a second order equation in t, requires the specification of the initial velocity.
  • The boundary conditions are insufficient.  For the biharmonic operator you need to specify not only the values of the unknown on the boundary, but also its derivatives.

Maple 2020 is unable to solve the following much simpler biharmonic equation in 1D:

restart;
pde := diff(w(x,t), t$2) + diff(w(x,t), x$4) = 0;
bc := w(0,t) = 0, w(1,t) = 0, D[1](w)(0,t)=0, D[1](w)(1,t)=0;
ic := w(x,0) = sin(Pi*x)^2, D[2](w)(x,0)=0;
pdsolve({pde,ic,bc});

I don't have Maple 2021. Give it a try with this 1D problem first.

PS: Maple 2020 can easily solve the 1D  problem above numerically. That doesn't help you with your 2D problem since Maple 2020's numeric PDE solver cannot handle PDEs with three independent variables.

 

 

The antiderivative of a function is determined up to an arbitrary additive constant.  So saying you want the antiderivative to be 0.739 is not meaningful.

Here is what you want:

restart;
f := sqrt(4+1/x)/2;
integ := int(f, x=0..a) assuming a > 0;
ans := fsolve(integ=0.739);
                      ans := 0.3726369981
# verify the answer
int(f, x=0..ans);
                          0.7390000000

 

As to your question about the shape of a chain draped over two pulleys, I suggest that you look first at the much simpler problem of determining the shape of a chain of length L suspended from two points within a vertical plane. Here is how it is done.

Let (x1,y1) and (x2,y2) be the Cartesian coordinates of the suspension points, and let y = f(x) be the equation of the chain.

The curve goes through the points (x1,y1) and (x2,y2).  That gives us two equations:
eq1 := y1 = f(x1),
eq2 := y2 = f(x2).

The curve's length is calculated by the usual curve length formula from calculus:
eq3 := int( sqrt(1 + f'(x)^2), x=x1..x2) = L;

Okay, so far we have specified the chain's geometry.  Now the physics.  It's possible to show, but I will not go through it here, that the chain's static equilibrium is expressed through the differential equation f''(x) = k*sqrt(1+f'(x)^2), where k is to be determined. In the worksheet below we see that this differential equation's general solution is a hyperbolic cosine.

Solving the differential equation brings along two integration constants c1 and c2.  So altogether, the differential equation's solution involves three unknown c1, c2, k.  We solve the system of equations {eq1, eq2, eq3} to determine those three unknowns.

The system is transcendental and has no useful symbolic solution, but we may solve the system numerically for any desired set of parameters (x1,y1), (x2,y2), and L.  See the worksheet below.

Going back to your question on the shape of a chain draped over two pulleys, that can be done in a similar way. The curve will be a catenary since that's the solution of the differential equation of the chain's equilibrium. I haven't done the calculations myself but the idea would be that instead of passing the curve from the points (x1,y1), (x2,y2), we need to find the curve which is tangent to the circles c1 and c2 that represent the pulleys. The solution should be within reach but not worth the trouble unless you really need it for some purpose.

Finally, regarding your question whether the points labeled A in your diagram in your worksheet are on horizontal diameters, the answer is no.  Just consider a chain which is tightly wrapped around the pulleys. The points A will be on vertical diameters in that case.

restart;
# The differential equation of a hanging chain
de := diff(y(x),x,x) = k*sqrt(1 + diff(y(x),x)^2);

# Solving the differential equation we obtain the equation of the hanging
# chain which involves three parameters, _C1, _C2, k:
dsolve(de):
f := unapply(rhs(%), x);

# We want the chain to go through the points (x1,y1) and (x2,y2):
eq1 := y1 = f(x1);
eq2 := y2 = f(x2);

# Additionally, we want the chain's length be a specified amount, L.
# The chain's length is calculated by the curve length formula from calculus:
int(sqrt(1+D(f)(x)^2), x):
simplify(%) assuming real:
eq3 := eval(%, x=x2) - eval(%, x=x1) = L;

# Now we need to solve the system of three equations {eq1, eq2, eq3} for the
# three unknowns _C1, _C2, k.  There is no symbolic solution, so we do  a
# numerical illustration.
params := { x1=0, y1=0, x2=10, y2=2, L=15 };
sys := eval({eq1,eq2,eq3}, params);
sol := fsolve(sys, {_C1,_C2,k}, {k=0..infinity});
#
# And now plot the solution
%plot(f(x), x=x1..x2, scaling=constrained):
eval(%, sol union params):
value(%);

 

Download worksheet: mw.mw

 

Here is the solution of your revised problem.  The solutions for f(x) and g(x) can be obtained entirely within Maple.

restart;

# Write the f derivative with the D operator
eq1 := x*D(f)(x^2) + diff(g(x),x) = cos(x) - 3*x^3;
eq2 := f(x^2) = sin(x) - x^4 - g(x);

# Calculate f'(x^2) from eq2 and substitute the result in eq1
# to obtain a differential equation in the unknown g(x):
diff(eq2, x);
isolate(%, D(f)(x^2));
de := subs(%, eq1);

# Solve the de for g(x), and then plug it in eq2 to calculate f(x);
dsolve(de);
eval(eq2, %);
algsubs(x^2=x, %);

The calculated f(x) and g(x) agree with what you have anticipated.

 

In your worksheet you refer to Problem #104 of 200 Puzzling Physics Problems. It is shown there that the shape of a chain in a steady motion between pulleys can be anything, so asking how to describe the chain's shape is pointless.

I must point out that the book's solution ignores gravity.  Adding gravity will drastically change the problem and make it quite nontrivial.  I don't expect a simple/elegant solution in that case.

The cantenary solution noted by Carl is good for a static chain.  As he has noted, a static uniform chain supported by two uprights of equal height assumes the shape of a simple hyperbolic cosine function. In fact, the "equal heights" assumption is unnecessary; the shape is a hyperbolic cosine even if the supports are at different levels.

 

Your candidates for f(x) and g(x) will work if you change x*f'(x) to x^2*f'(x) in the first equation:

restart;
eq1 := x^2*diff(f(x),x) + diff(g(x),x) = cos(x) - 3*x^3;
eq2 := f(x^2) + g(x) = sin(x) - x^4;
f := x -> -1/2*x^2 + C;
g := x -> sin(x) - 1/2*x^4 - C;

(lhs-rhs)(eq1);
                               0
(lhs-rhs)(eq2);
                               0

The above verifies that your candidates for f and g are correct.  I don't know how to solve the equations if f and g were not given. I don't expect that to be doable in Maple.

First 15 16 17 18 19 20 21 Last Page 17 of 58