Rouben Rostamian

MaplePrimes Activity


These are replies submitted by Rouben Rostamian

As Pólya is reported to have said: For any problem that you cannot solve, there is an easier problem that you cannot solve.  Solve that one first!

OK, so here is an easier problem obtained by setting μ=0 and reversing the 0 and 1 in the initial condition:

restart;
pde := diff(u(x,t),t) + u(x,t)*(diff(u(x,t),x)) = 0;
ic := u(x,0) = piecewise(x <= 0, 0, 1);

This is the classic Riemann problem.  Its exact solution is u(x,t)=F(x/t), where
F := s -> piecewise(s<=0, 0, s<1, s, 1);

I am unable to coax Maple's pdsolve() to arrive at that solution.  Perhaps someone else can. I would have expected HINT=f(x/t) to do the job but it doesn't.

 

@devraj You need to supply an initial condition as well.  I will take c(y,0) = cos(Pi*y).  Change as needed.

Now you may solve the problem numerically as follows:

restart;
pde := Omega*(diff(c(y, t), t))-Omega*Lambda*sin(t) = diff(c(y, t), y, y);
bc := (D[1](c))(0, t) = 0, (D[1](c))(1, t) = 0;
ic := c(y,0) = cos(Pi*y);
params := Omega=100, Lambda=0.001, pe=10;
sol := pdsolve(eval(pde, {params}), {bc, ic}, numeric);
sol:-plot3d(y=0..1, t=0..30, style=patchcontour);

This problem is very different from that you originally asked, so it's doubtful whether this solution is of any value.

 

@devraj Let's look at your new set of equations (1) and (2).

Differentiating equation (1) with respect to x leaves us with du/dx=0.  That says the term du/dx in equation (2) is zero and should be dropped.

Somehow I don't think that this is what you are looking for.  Ask someone who knows more about this to help you.

@devraj You should have a closer look at your equations because as they are, they make no sense at all.

Let's take, for example, the equations that you have listed as (1.) and (2.) in your post.  Subtract them.  You get:
u = pe*diff(u(x,y), x) .

On the other hand, you also say:

u=6*Lambda*(-y^2+y)*cos(2*Pi*x)+1

These two expressions for u directly contradict each other.

I suggest that you forget about Maple for now.  Figure out what your equations should be, or better yet, talk with someone who can guide you with those.  Once you have a logically consistent problem, come back here and then we will see how Maple can help.

 

@devraj You forgot to answer my question "What is x?".  I am referring to the cos(2*Pi*x) that appears in your PDE.

@vv Ah, now that's better :-)

@vv That expression looks nice but unfortunately it's not very practical as can be seen when trying
plot(answer, x=Pi/2-0.2..Pi/2+0.2);

@nm  In this problem the initial and boundary conditions are artificially designed/fine-tuned to yield a simple analytical solution.  The slightest change, for instance replacing the 12 by 13, will destroy the balance and the simple solution will be lost.  In the general case the solution would be expressed as an infinite sum—a Fourier series.

As to your question on how I thought of the "+" form of the solution, I looked at the boundary conditions and saw that the two ends move together, one end moving as −12 t2, and the other lagging by 1 relative to the other.  That led me to think of a rigid translational motion.  Since the initial condition is given as x4, the natural candidate would be x4 − 12 t2.  I tried that in my head, and saw that it works.

That gave me the idea of the HINT=`+` in Maple, and that yielded the same result.

This does not make the HINT=`+` option a useful choice for such problems.  That choice works in this case by "accident" because of the artificially fine-tuned conditions.

@assma Did you try saving it as EPS?

@assma I don't have Maple 18 to test, but it is likely that it should be able to handle the saving of a 3D figure as EPS.

Note, however, that the EPS format cannot handle transparency.  If your figure has transparent parts, then you will have to save it in the JPG format which is not the best thing since it is a lossy format and results in somewhat blurry images.

@josephrajck The equation and the boundary conditions are correct, but Initial conditions are missing.  Additionally, you can't get a numerical solution if you leave the coefficient "a" and the length "l" inspecified.  Here I will take both of these to be 1.  You may change then as needed.

restart;

pde := diff(y(x,t),t,t) + diff(y(x,t),x$4) = 0;

diff(diff(y(x, t), t), t)+diff(diff(diff(diff(y(x, t), x), x), x), x) = 0

bc := y(0,t)=0,
      D[1](y)(0,t)=0,
      D[1,1](y)(1,t)=0,
      D[1,1,1](y)(1,t)=0;

y(0, t) = 0, (D[1](y))(0, t) = 0, (D[1, 1](y))(1, t) = 0, (D[1, 1, 1](y))(1, t) = 0

ic := y(x,0) = x^2,
      D[2](y)(x,0)=0;

y(x, 0) = x^2, (D[2](y))(x, 0) = 0

sol := pdsolve(pde, {bc, ic}, numeric);

_m140082965461792

sol:-animate(t=5.23, frames=100, thickness=12, title="time = %f");

Download mw.mw

@vv Uh, OK, I see it now and I agree.  Also thanks to Preben for his explanations.

 

@Preben Alsholm Aren't all variables at the top level global by default?
restart;
type(a, `global`);          #true
type(whatever, `global`);   #true
type(_Z, `global`);         #true

I don't understand what it is so special about singling out _Z as being global.

 

@Preben Alsholm The error message that you have shown does not complain about "global".  It complains about "protected".  It may be a global variable, but that's beside the point.

In fact _Z is not unique in that respect.  Try any of:

_a := 12;
_b := 12;
_c := 12;

In general it's good to work under the assumption that all identifiers beginning with underscore are reserved for Maple's internal use, and avoid messing around with them.  Even if some of such identifiers are unprotected now, they may become unavailable in future releases.

This is exactly parallel to the C programming language where identifiers beginning with underscore are reserved for the compiler's internal use.  One may assign to them at one's own risk.

 

I haven't checked all of the calculations in that manuscript, however just

spot-checking it reveals some errors which weakens my trust in the rest

of its statements.

 

First, I have already noted the F(t, x, y) versus F(t, x, y, z) confusion,

but perhaps that can be attributed to a careless copy-paste error.  It is

clear that all F(t, x, y, z) and G(t, x, y, z) are supposed to be F(t, x, y)

and "G(t,x,y)."

 

Second, each of the two phase portraits in Figure 1 on page 161 show

an equilibrium near x = .5, y = 0.6e-1.  That is inconsistent with the

given data as we see in the calculations below. Again, that may be due

to typographical errors in the article but there is no way for us to tell.

 

Of course, it is possible that I am mistaken, in which case I would like

to be corrected.

 

restart;

The system of differential equations is defined in equation (1.2) in terms

of the functions F and G:

F := r*x*(1-x/k) - beta*x*y/(a+x^2);
G := mu*beta*x*y/(a+x^2) - d*x - eta*x*y;

r*x*(1-x/k)-beta*x*y/(x^2+a)

 

mu*beta*x*y/(x^2+a)-d*x-eta*x*y

(1)

The numerical values of the coefficients are specified on page 159:

params := r=0.05, a=0.8, mu=0.8, d=0.24,
          eta=0.01, beta=0.6, k=1.6;

r = 0.5e-1, a = .8, mu = .8, d = .24, eta = 0.1e-1, beta = .6, k = 1.6

(2)

The equilibria are the roots of the system F = 0, G = 0:

sys := eval({F=0, G=0}, {params});

{0.5e-1*x*(1-.6250000000*x)-.6*x*y/(x^2+.8) = 0, .48*x*y/(x^2+.8)-.24*x-0.1e-1*x*y = 0}

(3)

Figure 1 on page 161 indicates that there is an equilibrium at around "(x=0. 5, y=0.06)".

However we may verify in a number of ways that there are no equilibria near

that point at all.  One (but not the only) way is through Maple's fsolve() which fails to

find a root in the specified region:

fsolve(sys, {x=0.1..1.0,y=0.01..0.10});

fsolve({0.5e-1*x*(1-.6250000000*x)-.6*x*y/(x^2+.8) = 0, .48*x*y/(x^2+.8)-.24*x-0.1e-1*x*y = 0}, {x, y}, {x = .1 .. 1.0, y = 0.1e-1 .. .10})

(4)
 

 

Download mw2.mw

First 61 62 63 64 65 66 67 Last Page 63 of 99