Advanced Search

Rouben Rostamian

8501 Reputation

23 Badges

19 years, 31 days
Contact Rouben Rostamian

MaplePrimes Activity


These are replies submitted by Rouben Rostamian

The solution without a source term...

Rouben Rostamian 8501
11 hours ago
0 0

For whatever it's worth, here is the solution of your PDE in the absense of a source terms. 
 

> 

restart;

Bipolar to Cartesian transformation F is defined as

> 

F := (tau,sigma) -> [
        a*sinh(tau)/(cosh(tau)-cos(sigma)),
                a*sin(sigma)/(cosh(tau)-cos(sigma)) ];

proc (tau, sigma) options operator, arrow; [a*sinh(tau)/(cosh(tau)-cos(sigma)), a*sin(sigma)/(cosh(tau)-cos(sigma))] end proc

where a is an arbitrary parameter, and tau = -infinity .. infinity, and sigma = -Pi .. Pi.NULL

The rectangle

Omega = `&x`({`&tau;__1` < tau and tau < `&tau;__2`}, {-Pi < sigma and sigma < Pi})

in the tau, sigma plane is mapped onto the region omegabetween two off-center circles
in the Cartesian plane.

Specifically, the coordinate lines "tau&equiv;"constant map into circles centered
on the x axis in the Cartesian plane, and in particular, the images of the
rectangle's edges tau = `&tau;__1` and tau = `&tau;__2`map to the bounding circles of omega.

 

Let's calculate the radius and center of the circle corresponding to an

arbitrary "tau."

> 

F(tau,  0):  # the x coordinates of where
F(tau, Pi):  # the circle intersects the x axis
simplify(1/2*(%% - %)[1]):
radius := unapply(%, tau);

proc (tau) options operator, arrow; a*csch(tau) end proc

> 

F(tau,  0):  # the x coordinates of where
F(tau, Pi):  # the circle intersects the x axis
simplify(1/2*(%% + %)[1]):
center := unapply(%, tau);

proc (tau) options operator, arrow; a*coth(tau) end proc

We wish to chose the three parameters a, `&tau;__1`, `&tau;__2` so that the inner and outer
boundaries of omega are circles of given radii r and R, and the center of

the inner circle is offset relative to the outer circle by the amount kappa*R.

That results in three equations which we solve for the unknowns

a, `&tau;__1`, `&tau;__2`.  Maple returns a sequence of five solutions, the first four are

not interesting.  We pick the last one.

> 

radius(tau__2) = R,
radius(tau__1) = r,
center(tau__2) - center(tau__1) = kappa*R:
solve({%}, {a, tau__1, tau__2}):
params_general := allvalues(%[-1])[1];

{a = (1/2)*((1/2)*R^2*kappa^3*(R^2*kappa^2+R^2-r^2+(R^4*kappa^4-2*R^4*kappa^2-2*R^2*kappa^2*r^2+R^4-2*R^2*r^2+r^4)^(1/2))-R^4*kappa^3-(1/2)*R^2*kappa*(R^2*kappa^2+R^2-r^2+(R^4*kappa^4-2*R^4*kappa^2-2*R^2*kappa^2*r^2+R^4-2*R^2*r^2+r^4)^(1/2))-kappa*r^2*(R^2*kappa^2+R^2-r^2+(R^4*kappa^4-2*R^4*kappa^2-2*R^2*kappa^2*r^2+R^4-2*R^2*r^2+r^4)^(1/2))+R^4*kappa+R^2*kappa*r^2-(1/2)*r^2*(R^2*kappa^2+R^2-r^2+(R^4*kappa^4-2*R^4*kappa^2-2*R^2*kappa^2*r^2+R^4-2*R^2*r^2+r^4)^(1/2))/kappa+(1/2)*r^4*(R^2*kappa^2+R^2-r^2+(R^4*kappa^4-2*R^4*kappa^2-2*R^2*kappa^2*r^2+R^4-2*R^2*r^2+r^4)^(1/2))/(R^2*kappa))/(R*kappa*((1/2)*kappa*(R^2*kappa^2+R^2-r^2+(R^4*kappa^4-2*R^4*kappa^2-2*R^2*kappa^2*r^2+R^4-2*R^2*r^2+r^4)^(1/2))-R^2*kappa-(1/2)*r^2*(R^2*kappa^2+R^2-r^2+(R^4*kappa^4-2*R^4*kappa^2-2*R^2*kappa^2*r^2+R^4-2*R^2*r^2+r^4)^(1/2))/(R^2*kappa))), tau__1 = ln(-(1/2)*(R^2*kappa^2+R^2-r^2+(R^4*kappa^4-2*R^4*kappa^2-2*R^2*kappa^2*r^2+R^4-2*R^2*r^2+r^4)^(1/2))*r/(R^3*kappa*((1/2)*(R^2*kappa^2+R^2-r^2+(R^4*kappa^4-2*R^4*kappa^2-2*R^2*kappa^2*r^2+R^4-2*R^2*r^2+r^4)^(1/2))/R^2-1))), tau__2 = ln((1/2)*(R^2*kappa^2+R^2-r^2+(R^4*kappa^4-2*R^4*kappa^2-2*R^2*kappa^2*r^2+R^4-2*R^2*r^2+r^4)^(1/2))/(R^2*kappa))}

That's explicit but rather unpleasant.  Let's plug in specific numbers for R, r, kappa:

> 

params := simplify(rationalize(
        subs(R=1, r=1/2, kappa=1/4, params_general)));

{a = (1/8)*105^(1/2), tau__1 = -2*ln(2)+ln(11+105^(1/2)), tau__2 = -3*ln(2)+ln(105^(1/2)+13)}

Now that's not bad.  Let's plot the region to see what omega looks like;

> 

subs(params, F(tau__1, sigma)):
p1 := plot([%[], sigma=-Pi..Pi], scaling=constrained, color=red):
subs(params, F(tau__2, sigma)):
p2 := plot([%[], sigma=-Pi..Pi], scaling=constrained, color=blue):
plots:-display(p1,p2);

We wish to solve Laplace's equation in omega with the boundary conditions u = 0

on the outer circle and u = 1 on the inner circle.  But omega is the image of the

rectangle Omega, and Laplace's equation is invariant under the mapping, therefore

it suffices to solve the equation on the rectangle with boundary conditions

u = 0 on the edge tau = `&tau;__2` and u = 1 on the edge tau = `&tau;__1`, and periodic boundary

conditions on the other two edges.  That's straightforward to solve by hand

and we find that the solution is linear in tau:

tau

> 

U := 1/(tau__1 - tau__2)*(tau - tau__2);

(tau-tau__2)/(tau__1-tau__2)

Here is the graph of the temperature in Cartesian coordinates:

> 

%plot3d([F(tau,sigma)[], U], tau=tau__1..tau__2, sigma=-Pi..Pi,
        shading=zhue, style=patchcontour):
subs(params, %):
value(%);

 

Download heat-equation-in-off-center-annulus.mw

 

Graphs...

Rouben Rostamian 8501
July 22 2025
0 0

@Paras31 I see that your graphs of Short-time behavior are noticeably different from the book's.  That may be due to numerical errors that may have crept in.

The following graphs, plotted using Maple's solution that I had noted in my previous reply, seem to be identical to the book's.

Comments...

Rouben Rostamian 8501
July 21 2025
0 0

You have w0 and w[0].  Are these supposed to be the same thing?

You have w[0](eta,z).  What is w[0]?

You have w(r,z) := sum(p^i*w[i](r, z), i = 0 .. N); 
This is attempting to define w in terms of w.  Makes no sense.

You have D*w[0](0,z).  I am guessing that the D here is supposed to be the differentiation operator. But w[0] is a function of two variables, so you need to say derivative with respect to which variable.

Conclusion: What you are attempting to do is too complex for a first project in Maple.  You should spend some time learning the basics of the language before attempting this.

 

No Maple for Screen Readers for Linux...

Rouben Rostamian 8501
June 30 2025

@C_R I have not seen Maple for Screen Readers because it is specific to Windows, which I don't have.  I will give it a try if/when Maplesoft comes up with a Linux version.

Corrections...

Rouben Rostamian 8501
June 28 2025
0 0

@Math-dashti In my worksheet I had changed the symbol "eta" of your textbook's equation (7.1) to "t" because t is easier to type.  Your modified worksheet does not work because it now mixes up t and eta.

The next image of the textbook's page that you have posted, ends with "... can be simplified to the following nonlinear ODE".  The ODEs, which come in the next page, are not shown.  You need to replace the ODEs in my worksheet with those ODEs.

And as you do that, you need to think about what to do with "x" that enters the equations.  You cannot solve the equations and obtain a plot without knowing what x is.  I don't know the context, so I cannot help you there.

The complete phase portrait...

Rouben Rostamian 8501
June 26 2025
0 0

@Math-dashti This may be what you want.

> 

restart;
> 

with(plots):

Let's write diff(u(x), x) = f(u(x), v(x)), diff(v(x), x) = g(u(x), v(x)) for the two differential equations.

> 

f := (u,v) -> v;
g := (u,v) -> F__1*u + F__2*u^2;

proc (u, v) options operator, arrow; v end proc

proc (u, v) options operator, arrow; F__1*u+F__2*u^2 end proc

The equilibria:

> 

equilibria := solve({f(u,v)=0, g(u,v)=0}, {u,v});

{u = 0, v = 0}, {u = -F__1/F__2, v = 0}

The differential equations (not needed in this worksheet)

> 

de1 := diff(u(t),t) = f(u(t),v(t));
de2 := diff(v(t),t) = g(u(t),v(t));

diff(u(t), t) = v(t)

diff(v(t), t) = F__1*u(t)+F__2*u(t)^2

The conserved quantity

> 

P := 1/2*F__1*u^2 + 1/3*F__2*u^3 - 1/2*v^2;

(1/2)*F__1*u^2+(1/3)*F__2*u^3-(1/2)*v^2

Let's take F__1 = 1, F__2 = -1NULL

> 

F__1, F__2 := 1, -1;

1, -1

Here are the equilibra now

> 

equilibria;

{u = 0, v = 0}, {u = 1, v = 0}

Plot the phase portrait

> 

p1 := contourplot(P, u=-1.5..2, v=-1.5..1.5, scaling=constrained,
    colorscheme="DivergeRainbow", contours=[seq](x, x=-0.4..0.4,0.1)):

Mark the equilibria:

> 

p2 := pointplot([[0,0],[1,0]], symbol=solidcircle, symbolsize=15, color=black):

Add direction field

> 

p3 := fieldplot([f(u,v), g(u,v)], u=-1.5..2, v=-1.5..1.5,
        arrows=medium, fieldstrength=fixed(0.4), grid=[10,10], labels=["",""]):

Put all pieces together

> 

display(p1,p2,p3, axes=normal);

 

Download phase-portrait.mw

 

Good code...

Rouben Rostamian 8501
June 25 2025

@dharr Considering the frequent requests for Adomian polynomials in this forum, your LAD code is can be a great time-saver.  I will certainly refer the readers to it when the occasion comes up.

@dharr Okay, perhaps I am looking a...

Rouben Rostamian 8501
June 25 2025

@dharr Okay, perhaps I am looking at this through the wrong lens.  I had expected the series solution to provide an alternative to a numerical solution.  It looks like that's not the case.

There may be some information content in the series solution but I have a hard time getting my hands on it.  Suppose I am interested in the porous medium equation with the initial condition u(x,0)=exp(−x^2) for which there is no symbolic solution.  What sort of useful information do I get out of the LAD series?

@TechnicalSupport Does that work fo...

Rouben Rostamian 8501
June 25 2025

@TechnicalSupport Does that work for you on Linux?  It doesn't for me.

Solution at x=0...

Rouben Rostamian 8501
June 25 2025

@dharr Here is what we get when we evaluate the series solution in LAD_Barenblatt.mw at x=0 and plot the result as a function of time:

> 

subs(x=0, approx);

1-(75/4)*t+(28125/32)*t^2-(6328125/128)*t^3+(6169921875/2048)*t^4-(1573330078125/8192)*t^5+(825998291015625/65536)*t^6-(221249542236328125/262144)*t^7+(481217754364013671875/8388608)*t^8

> 

plot(%, t=0..1);

We see that the series solution rises dramatically over time, while it's supposed to decrease to zero.

The non-differentiability of the solution does not seem to be the issue.  We get the same kind of incorrect behavior out of LAD if we replace the initial condition with the very well-behaved u(x,0)=exp(−x^2).

The remark in the Kybernetes article about the "consitency of boundary conditions" is not relevant here since we are not solving a boundary value problem; there are no boundaries as the spatial dimension is the entire x axis.

Wait until the next update...

Rouben Rostamian 8501
June 24 2025

@WD0HHU The stop button does nothing in Maple 2025 and Maple 2025.1.  Until that issue is fixed, I am sticking with Maple 2024.

@lucaud Before executing reflect(di...

Rouben Rostamian 8501
June 24 2025
0 0

@lucaud Before executing reflect(display...), you need to load the plots and plottools packages, as in

with(plots):
with(plottools):

If that doesn't work, then upload your entire worksheet for further comments.  To upload the worksheet, "reply" to your own original post, and within the toolbar of the editor window that comes up, click on the big fat green arrow.

That looks interesting...

Rouben Rostamian 8501
June 23 2025

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

Rouben Rostamian 8501
June 17 2025
0 0

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

dsolve's solution...

Rouben Rostamian 8501
June 17 2025
0 0

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

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