Rouben Rostamian

MaplePrimes Activity


These are answers submitted by Rouben Rostamian

In the first line of the image that you have attached, we see the input is
f := whatever
but the corresponding output is
f = whatever

Note the disagreement between the ":=" and "=" in the input versus the output.  This means that you have first executed the input without a colon and obtained the output, and then you have inserted the colon in the input but have not executed it!  As a result, f has not been assigned a value and therefore it appears as a plain f in the final result.

Execute the worksheet again without changing anything.  You should get the expected result.

 

 

Instead of plotting and then shifting, you may want to plot the arrows directly in the desired location.

For instance, to plot arrows originating at the point [3,4], you can do:

restart;
plots:-arrow([3,4], [[1,0], [0,1]]);

Moreover, as we see here, a single call to arrow is sufficient to plot any number of arrows.

With the unknown function u(x,y), the Dirichlet condition  u(0,y) = f(y) is entered as

u(0,y) = f(y);

The Neumann condition ∂u/∂x(0,y) = g(y) is entered as 

D[1](u)(0,y) = g(y);

The Neumann condition ∂u/∂y(x,0) = h(x) is entered as 

D[2](u)(x,0) = h(x);

The notations D[1] and D[2] refer to differentiation with respect to the first and second variables of u(x,y).

 

The presence of the cosine term is more of an effect of the geometry than physics, although I dislike putting a hard barrier between the two subjects.

The Pythagorean Theorem says that the square of the hypotenuse is the sum of the squares of the two other sides in a right triangle.  The generalization to an arbitrary triangle is known as the law of cosines.  I trust that you already know that but in the unlikely case, look it up in Wikipedia.

The law of cosines enters your calculations because the two vectors that you have shown are not perpendicular. The pearl's velocity is the sum of those two vectors.  Calculating the sum relies on the law of cosines.

restart;

Typesetting:-Settings(typesetdot):

 

# the horizontal vector labled a*(diff(phi(t), t))in the diagram
< a*diff(phi(t),t), 0 >:

# the slanted vector labeled a*(diff(theta(t), t)) in the diagram
a*diff(theta(t),t) * < -cos(theta(t)), sin(theta(t)) >:

# the sum of the two vectors above, which is the pearl's absolute velocity
v := %% + %;

Vector[column](%id = 18446884233678926542)

The pearl's kinetic energy

simplify(1/2*m*v^+ . v);

(1/2)*m*a^2*((diff(theta(t), t))^2-2*(diff(phi(t), t))*(diff(theta(t), t))*cos(theta(t))+(diff(phi(t), t))^2)

Download mw.mw

This works in Maple 2019.  It may require minor adjustments in older versions.

restart;

A := < 1, 2; 3, 4 >;

Matrix(2, 2, {(1, 1) = 1, (1, 2) = 2, (2, 1) = 3, (2, 2) = 4})

X := < x[1](t), x[2](t) >;

Vector(2, {(1) = x[1](t), (2) = x[2](t)})

sys := diff(X,t) =~ A . X;

Vector(2, {(1) = diff(x[1](t), t) = x[1](t)+2*x[2](t), (2) = diff(x[2](t), t) = 3*x[1](t)+4*x[2](t)})

dsolve(sys);

{x[1](t) = _C1*exp((1/2)*(5+33^(1/2))*t)+_C2*exp(-(1/2)*(-5+33^(1/2))*t), x[2](t) = -(1/4)*_C2*exp(-(1/2)*(-5+33^(1/2))*t)*33^(1/2)+(1/4)*_C1*exp((1/2)*(5+33^(1/2))*t)*33^(1/2)+(3/4)*_C2*exp(-(1/2)*(-5+33^(1/2))*t)+(3/4)*_C1*exp((1/2)*(5+33^(1/2))*t)}


 

Download mw.mw

As tomleslie said, there are many ways of achieving the result.  Here is a straightforwad "bare hands" way which may be more accessible to a beginner.

restart;

Angles of the two vectors with respect to the vertical, in radians:

a1 := convert(30.0*degrees, radians);
a2 := convert(60.0*degrees, radians);

.5235987757

1.047197551

The two vectors:

v1 := 1.0 * < sin(a1), cos(a1) >;
v2 := 1.5 * < sin(a2), cos(a2) >;

Vector(2, {(1) = .5000000001, (2) = .8660254037})

Vector[column](%id = 18446884770620611990)

Their sum:

v := v1 + v2;

Vector(2, {(1) = 1.79903810565, (2) = 1.6160254040000002})

v's length

sqrt(v[1]^2 + v[2]^2);

HFloat(2.418279597555689)

v's angle with respect to the vertical in radians and in degrees

arctan(v[1],v[2]);
convert(%, degrees);

HFloat(0.838936788996444)

HFloat(48.06753728147721)*degrees

 

 

 

Download mw.mw

Probably this is what you want:

for n in (10*2^i $i=0..4) do
    print(n);  # or whatever
end do;

or equivalently

for i from 0 to 4 do
    n := 10*2^i;
    print(n); # or whatever
end do;

or

n := 10; 
while n <= 160 do
    print(n); # or whatever
    n := 2*n;
end do;

Pick whichever you like.

 

 

Here is a possible way.

 

Example 1

 

restart;
ode:=3*diff(y(x),x)^2+diff(y(x),x)^3+sin(x)+y(x)=x*y(x)+x*diff(y(x),x);

3*(diff(y(x), x))^2+(diff(y(x), x))^3+sin(x)+y(x) = x*y(x)+x*(diff(y(x), x))

indets['flat'](ode,{`^`('identical'(diff(y(x),x)),'algebraic'),'identical'(diff(y(x),x))})

{(diff(y(x), x))^2, (diff(y(x), x))^3, diff(y(x), x)}

ode_wanted:= 3*diff(y(x),x)^2+diff(y(x),x)^3-x*diff(y(x),x)=-sin(x)-y(x)+x*y(x)

3*(diff(y(x), x))^2+(diff(y(x), x))^3-x*(diff(y(x), x)) = -sin(x)-y(x)+x*y(x)

selectremove(has, (lhs-rhs)~(ode), diff(y(x),x));
ode_new := %[1] = - %[2];
ode_new - ode_wanted;

3*(diff(y(x), x))^2+(diff(y(x), x))^3-x*(diff(y(x), x)), sin(x)+y(x)-x*y(x)

3*(diff(y(x), x))^2+(diff(y(x), x))^3-x*(diff(y(x), x)) = -sin(x)-y(x)+x*y(x)

0 = 0

 

 

Example 2

 

restart;
ode:=3*diff(y(x),x)^2+diff(y(x),x)=x*diff(y(x),x)+5;

3*(diff(y(x), x))^2+diff(y(x), x) = x*(diff(y(x), x))+5

indets['flat'](ode,{`^`('identical'(diff(y(x),x)),'algebraic'),'identical'(diff(y(x),x))})

{(diff(y(x), x))^2, diff(y(x), x)}

ode_wanted:= 3*diff(y(x),x)^2+diff(y(x),x)-x*diff(y(x),x)=5

3*(diff(y(x), x))^2+diff(y(x), x)-x*(diff(y(x), x)) = 5

selectremove(has, (lhs-rhs)~(ode), diff(y(x),x));
ode_new := %[1] = - %[2];
ode_new - ode_wanted;

3*(diff(y(x), x))^2+diff(y(x), x)-x*(diff(y(x), x)), -5

3*(diff(y(x), x))^2+diff(y(x), x)-x*(diff(y(x), x)) = 5

0 = 0

 


 

Download lhs_rhs-new.mw

 

I have made several changes to your worksheet and I hope that I have not introduced any errors.  Check to be sure.

The optimal diffusion coefficient turns out to be d = 0.14 (approximately).

restart;

t_number:=<0, 10, 20, 30, 40>:
m_number:=<11.50000000, 4.641182547, 1.273311993, 0.361845238, 0.288711649>:

 

Q:=proc(d)
    local pds, i, S, PDE, IBC,
          L := 2, Mx0 := 0.05, cx0 := Mx0/(1-Mx0),
          Mdb_i := m_number[1], ct0 := Mdb_i;

    if not type(d, numeric) then return 'procname(_passed)' end if;

    PDE := diff(C(x,t),t)=d*diff(C(x,t),x,x);
    IBC := {C(x,0)=ct0, C(0,t)=cx0, D[1](C)(L,t)=0};

    pds := pdsolve(PDE,IBC,numeric);
    S := 0;
    for i from 1 to 5 do
        # solution at the desired time

        pds:-value(t=t_number[i], output=listprocedure);
        eval(C(x,t), %);              # extract the C(x,t) at that time
        int(%, 0..2, numeric) / L;    # compute the average of C(x,t) at that time
        S += (% - m_number[i])^2;     # accumulate the residuals
     end do;  
     return sqrt(S);
end proc:

# Try out the proc

Q(0.1), Q(0.2), Q(0.5), Q(1.0);

HFloat(2.0125989798253805), HFloat(1.9268252656215228), HFloat(4.344399869795477), HFloat(4.74473396546933)

plot(Q(d), d=0.1 .. 0.3, numpoints=5, view=0..2);

# We conclude that the best choice of d is approximately 0.14

 


 

Download zz.mw

 

To plot the derivatives, add the output=operator option to your dsolve(), as in
p := dsolve(%, vars, numeric, output = operator);

Do odeplot() as before to plot the solutions.  To plot the derivative of y(t) (also those of x[1](t) and x[2](t)) do
eval(diff(y(t), t), Sys);
eval(%, p);
plot(%, t = 0 .. 1500);

As to the higher order derivatives, y'', y''', etc., these are not quite well-defined because the right-hand sides of your differential equations contain discontinuities, which means that y' is discontinuous.  You don't want to differentiate a discontinuous function, do you?

You can actually see the discontinuity in y' by changing the plot(%, t = 0 .. 1500) in what I have shown to plot(%, t = 500 .. 1500).

As to the plotting of the bar graph, I don't think it's difficult but I don't have the time to go into it right now.  Perhaps someone else will.

 

Of the two solutions returned by dsolve() one corresponds to positive y and the other to negative y.

restart;
de := diff(y(x),x) = abs(y(x)) + 1;
                         d                   
                  de := --- y(x) = |y(x)| + 1
                         dx                  
dsol := dsolve(de);
                       exp(-x)                           
      dsol := y(x) = - ------- + 1, y(x) = exp(x) _C1 - 1
                         _C1                             
odetest(dsol[1], de) assuming y(x)<0;
                               0
odetest(dsol[2], de) assuming y(x)>0;
                               0

Your differential equation has y(x) = 0 as a solution, and that's what
Maple returns by default.  However, for special choices of omega^2 there are
nonzero solutions.  Those special choices of omega^2 are called the problem's
eigenvalues, and their corresponding solutions are called the eigenfunctions.

Here is how we go about finding the eigenvalues and eigenfunctions.

We begin with a couple of self-evident observations.
   

1. 

If a function y(x) is an eigenfunction, then for any nonzero constant c
the function c*y(x) is also an eigenfunction.

2. 

The derivative "y '(0)" of an eigenfunction cannot be zero because the
conditions y(0) = 0 and "y '(0)=0" will imply that `&equiv;`(y(x), 0) which
is not possible since an eigenfunction is a nonzero function by definition.


Putting those two observations together, we may take "y '(0)=1" for
all eigenfunctions without a loss of generality.

This adds an extra boundary condition over and above what you
already have.  That gives us the flexibility of introducing a new
unknown in the system.  We take the new unknown to be the
constantomega and we solve the system for the unknowns {omega, y(x)}.

restart;

The differential equation

de := T*diff(y(x), x, x) + rho*omega^2*y(x) = 0;

T*(diff(diff(y(x), x), x))+rho*omega^2*y(x) = 0

The boundary conditions (note the extra third condition!)

bc := y(0) = 0, y(L) = 0, D(y)(0)=1;

y(0) = 0, y(L) = 0, (D(y))(0) = 1

dsol_tmp := dsolve({de,bc}, {y(x),omega});

{omega = Pi*(2*_Z1+_B1)*T^(1/2)/(rho^(1/2)*L), y(x) = L*sin(Pi*(2*_Z1+_B1)*x/L)/(Pi*(2*_Z1+_B1))}

Here omega and y(x) are expressed in terms of the arbitrary
constants _Z1 and _B1.  Let's see what they are:

about(_Z1);
about(_B1);

Originally _Z1, renamed _Z1~:

  is assumed to be: integer

Originally _B1, renamed _B1~:
  is assumed to be: OrProp(0,1)

OK then.  _Z1 is an integer and _B1 is either zero or one.  It follows
that the combination 2*_Z1+_B1 which enters the solution is an
arbitrary integer.  We call it n

dsol := eval(dsol_tmp, {2*_Z1+_B1=n});

{omega = Pi*n*T^(1/2)/(rho^(1/2)*L), y(x) = L*sin(Pi*n*x/L)/(Pi*n)}

Let's verify that this satisfies the differential equation

eval(de, dsol);

0 = 0


 

Download mw.mw

I don't understand what you mean by "returns only one option".  Here is the full solution.

restart;

wW := unapply(piecewise(0 <= x and x < L/2, 0, L/2 <= x and x <= L, w[0]), x);

wW := proc (x) options operator, arrow; piecewise(0 <= x and x < (1/2)*L, 0, (1/2)*L <= x and x <= L, w[0]) end proc

eq := k*diff(y(x), x$4) = wW(x);

eq := k*(diff(y(x), x, x, x, x)) = piecewise(0 <= x and x < (1/2)*L, 0, (1/2)*L <= x and x <= L, w[0])

dsolve({eq, y(0) = 0, D(y)(0) = 0, (D@@2)(y)(L) = 0, (D@@3)(y)(L)=0}, y(x))
  assuming 0 < L;

y(x) = piecewise(x < (1/2)*L, -L*w[0]*x^3/(12*k)+3*L^2*w[0]*x^2/(16*k), x < L, -L*w[0]*x^3/(6*k)+L^2*w[0]*x^2/(4*k)-L^3*w[0]*x/(48*k)+L^4*w[0]/(384*k)+w[0]*x^4/(24*k), L <= x, 7*L^3*w[0]*x/(48*k)-5*L^4*w[0]/(128*k))


 

Download mw.mw

In 2.mw we have
read "1.mw"; 
but Maple's read command is meant for reading plain text files, not Maple worksheets.

To fix, copy and paste the contents of 1.mw as a plain text into a file, let's say 1.maple.  Then change the read command to 
read "1.maple"; 

Here is 1.maple for your convenience.

 

 

Approximate solution

 

On dividing through by r, your differential equation changes to
d/(r*dr)*(r*d*H(r)/dr)+(k^2-b^2*r/R^2)*H(r) = 0.``
The first term is the second derivative of H(r) in polar coordinates.

With your numerical parameters, which are
R = 370, k = 100, b = 35,
the coefficient of the second term is 10000-(49/5476)*r, which varies
between 9996 and 10000 on the range of interest 1/R < r and r < R, and
therefore is essentially a constant. Consequently, for practical purposes
we may take b = 0 and reduce the differential equation to
d/(r*dr)*(r*d*H(r)/dr)+k^2*H(r) = 0,
We expect the solution to be something like cos(k*r), that is, cos(100*r),
and therefore highly oscillatory.  We conclude that searching for an
approximation in the form of power series is hopeless.

The good news is that Maple can obtain a symbolic solution for the
reduced differential equation.  Here it is.

restart;

de := 1/r*Diff(r*Diff(H(r),r),r) + k^2*H(r);

(Diff(r*(Diff(H(r), r)), r))/r+k^2*H(r)

bc := H(R)=0, D(H)(1/R)=R;

H(R) = 0, (D(H))(1/R) = R

dsol := dsolve({de,bc}, H(r));

H(r) = -R*BesselY(0, k*R)*BesselJ(0, k*r)/(k*(BesselJ(1, k/R)*BesselY(0, k*R)-BesselY(1, k/R)*BesselJ(0, k*R)))+BesselJ(0, k*R)*R*BesselY(0, k*r)/(k*(BesselJ(1, k/R)*BesselY(0, k*R)-BesselY(1, k/R)*BesselJ(0, k*R)))

Apply the parameters and plot:

R := 370:  k := 100:

plot(rhs(dsol), r=1/R..R, view=-0.1 .. 0.1, numpoints=10000);

We cannot see the graph's details because it oscillates so much, but we do see
the details if we plot it on a narrow range:

plot(rhs(dsol), r=100 .. 101);

which shows fast oscillations as expected.


 

Download de2.mw

First 21 22 23 24 25 26 27 Last Page 23 of 53