Rouben Rostamian

MaplePrimes Activity


These are replies submitted by Rouben Rostamian

I have occasionally heard about Adomian's method but have not had the occasion to delve into it in any depth, and to be honest, I have harbored some doubts about the utility of it.

But I will definitely be a convert if it can produce anything resembling the correct solution in the case of the classical porous medium equation.  I tried your proc but did not get anything useful, although it is very well possible that I did not apply it correctly.

Here is the porous medium equation along with an exact solution.  You may be able to do somehing with it.

PS: Perhaps the method will have a better chance of success with the initial condition u(x,0)=exp(-x^2).  There is no symbolic solution then, but pdsolve,numeric shows that the solution qiuckly takes the shape of an expanding semi-ellipse similar to the one in the animation below.

Barenblatt's solution to the porous medium equation

restart;

The porous medium equation with cubic nonlinearity

pde := Diff(u(x,t),t) = Diff(u(x,t)^3, x ,x);

Diff(u(x, t), t) = Diff(u(x, t)^3, x, x)

An initial condition in the form of the upper half of an ellipse:

ic := u(x,0) = 75^(1/4)*sqrt(max(0, sqrt(3)/15 - x^2*sqrt(75)/12));

u(x, 0) = 75^(1/4)*max(0, (1/15)*3^(1/2)-(5/12)*x^2*3^(1/2))^(1/2)

Exact solution, due to Barenblatt (1952):

sol := u(x,t) = sqrt(max(0, sqrt(3)/15 - x^2/(12*sqrt(t + 1/75))))/(t + 1/75)^(1/4);

u(x, t) = max(0, (1/15)*3^(1/2)-(5/4)*x^2/(225*t+3)^(1/2))^(1/2)/(t+1/75)^(1/4)

At any time t, the graph of the solution is the upper half of an ellipse.

The area under the ellipse is a constant Pi/5 at all times.

plots:-animate(plot, [rhs(sol), x=-2..2, color=red, thickness=4], t=0..4);

 

Download pme.mw

 

@janhardo The code that you have posted is good.  It implements what I called an alternative approach in my yesterday's message.

@vv Before seeing Carl Love's solution, I would have completely agreed with you.  I am impressed by seeing dsolve's solution pointed out by Carl.

@mmcdara Thanks for your analysis which shows how to obtain an approximate solution by smoothing out the transition at the interface.  In the absence of an explicit symbolic solution, that would have been a good alternative.

@Carl Love That's excellent, and works beyond my expectations.  Thanks!

@vv The solution that I have provided is called a weak solution of the boundary value problem.  The definition is technical so I would not go through it here, but they are the main objects of study in the functional analytical approach to PDEs.  In the case of the ODE which is the subject of this post, the weak solution ensures that the temperature, u(x), and the flux, a(x)u'(x), are continuous across the interface at x=0.  The derivative u'(x) is understood in the sense of distributions, or more precisely, in the sense of Sobolev spaces. The continuity of the flux amounts to the statement of conservation of thermal energy at x=0.

An alternative "by hand" solution may be obtained by finding the general solutions of the ODE separately on the x<0 and x>0 regions, and then determining the resulting four arbitrary constants by applying the boundary conditions at x=−1 and x=+1, as well as the continuity of the temperature and flux at x=0.

@John May Thanks for looking into this and offering a workaround.

@Math-dashti There is some flexibility in selecting the initial conditions, but they are not quite arbitrary.

Begin with finding the det, tra, dsc, and eigenvalues as I have shown in the worksheet.  With the help of that information, and the table and Figure 2.3.7 that you have provided, determine which of the ten possible phase portraits corresponds to your case.

Let's say that you have determined that the phase portrait is the one at the top-left of Figure 2.3.7. Looking at that phase portrait, you should be able to tell which initial conditions to pick.  That's how I picked the initial conditions in my solution of Exercise 7.

I suggest that you apply this method to solve Exercises 6, 4, 1 (in that order!) before attempting the others.

@C_R Maple's plottools:-rotate makes it possible to rotate a 3D plot about a vector n by angle phi if the objective is merely to rotate a PLOT3D structure.  To find the equation of the rotated torus, however, as it was done in the example, it will help to have something like quaternion algebra.  Furthermore, we will need full-fledged support for quaternion algebra and calculus to formulate differential equations of motion of a rigid body.  I will be happy to see that become reality in Maple some day.

@rlewis  It's difficult to provide help in the absence of concrete data.  As dharr pointed out, you can help others to help you if you upload your maple worksheet.

In the browser, hit Reply to your original post.  In the editor window that comes up, click on the big fat green arrow and follow instructions to upload your worksheet.

@acer Thanks for confirming the issue.  It did test this on Maple 2024 before posting, but did not go further back.  I will file an SCR.

@dharr Thanks for the workaround.  That's very helpful.

A lesson without an example is not very useful.  Here is a concrete example that illustrates how to rotate an object via quaternions.

restart;

with(plots):

The  product of the quaternions p=[a,u] and q=[b,v]

qProd := proc(p, q)
        local a, b, u, v;
        a := p[1];
        u := p[2];
        b := q[1];
        v := q[2];
        return [a*b - u^+.v, a*v + b*u + LinearAlgebra:-CrossProduct(u,v)];
end proc:

The conjugate of the quaternion p

qConj := proc(p)
        return [p[1], -p[2]];
end proc:

Vector(3, {(1) = 1/3, (2) = 2/3, (3) = 2/3})

(1/3)*Pi

Quaternion for rotating about the unit vector n by angle phi

n := < 1,2,2 >/3;
phi := Pi/3;
q := [ cos(phi/2), n*sin(phi/2) ];

Vector(3, {(1) = 1/3, (2) = 2/3, (3) = 2/3})

phi := (1/3)*Pi

[(1/2)*3^(1/2), Vector[column](%id = 36893622606788529436)]

Rotate the point r via the formula q o [0,r] o q*;

r := < x, y, z >;
rot := simplify(qProd(q, qProd([0,r], qConj(q))));

Vector(3, {(1) = x, (2) = y, (3) = z})

[0, Vector[column](%id = 36893622606788518484)]

The scalar part of the result is zero, as stated before.   The vector part of the result is the rotated position of r, let's call it v:

v := rot[2];

Vector(3, {(1) = (-(1/3)*y+(1/3)*z)*sqrt(3)+(5/9)*x+(1/9)*y+(1/9)*z, (2) = ((1/3)*x-(1/6)*z)*sqrt(3)+(1/9)*x+(13/18)*y+(2/9)*z, (3) = (-(1/3)*x+(1/6)*y)*sqrt(3)+(1/9)*x+(2/9)*y+(13/18)*z})

We see that v is related to r through a matrix multiplication, that is, v = R . r.

Here is the matrix R:

R := LinearAlgebra:-GenerateMatrix([seq](v), [x,y,z])[1];

Matrix(3, 3, {(1, 1) = 5/9, (1, 2) = -(1/3)*sqrt(3)+1/9, (1, 3) = (1/3)*sqrt(3)+1/9, (2, 1) = (1/3)*sqrt(3)+1/9, (2, 2) = 13/18, (2, 3) = -(1/6)*sqrt(3)+2/9, (3, 1) = -(1/3)*sqrt(3)+1/9, (3, 2) = (1/6)*sqrt(3)+2/9, (3, 3) = 13/18})

As a concrete example, consider the parametrically defined torus

T := < (3 + cos(t))*cos(s), (3 + cos(t))*sin(s), sin(t) >;

Vector(3, {(1) = (3+cos(t))*cos(s), (2) = (3+cos(t))*sin(s), (3) = sin(t)})

Rotate that torus about the vector n by angle phi:

T__rot := R . T;

Vector(3, {(1) = (5/3+(5/9)*cos(t))*cos(s)+(-(1/3)*sqrt(3)+1/9)*(3+cos(t))*sin(s)+((1/3)*sqrt(3)+1/9)*sin(t), (2) = ((1/3)*sqrt(3)+1/9)*(3+cos(t))*cos(s)+(13/6+(13/18)*cos(t))*sin(s)+(-(1/6)*sqrt(3)+2/9)*sin(t), (3) = (-(1/3)*sqrt(3)+1/9)*(3+cos(t))*cos(s)+((1/6)*sqrt(3)+2/9)*(3+cos(t))*sin(s)+(13/18)*sin(t)})

plot3d([T, T__rot], s=-Pi..Pi, t=-Pi..Pi, color=["Green","Orange"], style=surface, scaling=constrained);

 

Download rotation-quaternion.mw

That's rather odd.  If the point are collinear, then the area would be simply zero.  We don't need conditions for that.  Maple's answer would have made sense before the invention of zero.

@one man Okay, now I see what you are doing.  You specify the paths of the contact point independently on the two surfaces, and then roll one surface against the other while requiring the contact point to stay on both paths.  In that case there is one degree of freedom, and all is well.

As to the Latinization of Dragilev, I see that "Dragilev's Method" is written "Метод Драгилева" in Russian, so the Latinized name is "Dragilev".  There is no reason to insert an "h" there.

Edit: The impresario of Ballets Russes, Sergei Diaghilev (Сергей Дягилев), spelled his Latinized name with a "gh", but that was influenced by the Italian orthography, as he lived in Italy, died in Venice, and is buried there.

1 2 3 4 5 6 7 Last Page 2 of 99