Rouben Rostamian

MaplePrimes Activity


These are replies submitted by Rouben Rostamian

@denkasyan Your expression is different from what I have given.  Don't retype; just copy and paste my code into your worksheet.

@denkasyan Perhaps this is what you are asking:

restart;
with(PolyhedralSets):
y := -20;
p := PolyhedralSet([-x >= 0, -y >= 0, -z >= 0, -(1/2)*x >= 0, (-x-y)*(1/2) >= 0, (-x-z)*(1/2) >= 0]);
Plot(p);

If that's not what you want, then you will have to explain more clearly what it is that you are after.

 

@denkasyan  Saying "y = -20 is a constant" does not tell me what it is meant to do. How is it related to the polyhedron?

@Math Pi Euler Mariusz Iwaniuk has already shown one way of obtaining the higher order terms. Here I provide a more straighforward alternative method.

Below, you may change the value of N to obtain higher terms of any order.

restart;

de := diff(z(x),x,x) + z(x) - cos(2*x)/(1+delta*z(x)) = 0;

diff(diff(z(x), x), x)+z(x)-cos(2*x)/(1+delta*z(x)) = 0

 bc := z(-Pi/4)=0, z(Pi/4)=0

z(-(1/4)*Pi) = 0, z((1/4)*Pi) = 0

We are going to expand the solution of this boundary value problem in powers of delta,

and pick the leading N terms like this:

N := 3;  # change as needed
S := z(x) = add(u[n](x)*delta^n, n=0..N-1);

3

z(x) = u[0](x)+u[1](x)*delta+u[2](x)*delta^2

Substitute the series expression S in the differential equation and simplify:

eval(de, S):
series(lhs(%), delta, N):
collect(%, delata, simplify):
tmp := convert(%, polynom);

diff(diff(u[0](x), x), x)+u[0](x)-cos(2*x)+(diff(diff(u[1](x), x), x)+u[1](x)+cos(2*x)*u[0](x))*delta+(diff(diff(u[2](x), x), x)+(-u[0](x)^2+u[1](x))*cos(2*x)+u[2](x))*delta^2

This is supposed to be zero for all delta, therefore the coefficient of each power of
delta is zero.  That gives us a system of N coupled differential equations in N unknowns:

coeffs(tmp, delta);

diff(diff(u[0](x), x), x)+u[0](x)-cos(2*x), diff(diff(u[2](x), x), x)+(-u[0](x)^2+u[1](x))*cos(2*x)+u[2](x), diff(diff(u[1](x), x), x)+u[1](x)+cos(2*x)*u[0](x)

We extract boundary conditions for each u__n from the originally supplied boundary conditions,

and associate them with each of the differential equations obtained above:

sys := seq({coeff(tmp, delta, n), subs(z=u[n], [bc])[]}, n=0..N-1);

{diff(diff(u[0](x), x), x)+u[0](x)-cos(2*x), u[0](-(1/4)*Pi) = 0, u[0]((1/4)*Pi) = 0}, {diff(diff(u[1](x), x), x)+u[1](x)+cos(2*x)*u[0](x), u[1](-(1/4)*Pi) = 0, u[1]((1/4)*Pi) = 0}, {diff(diff(u[2](x), x), x)+(-u[0](x)^2+u[1](x))*cos(2*x)+u[2](x), u[2](-(1/4)*Pi) = 0, u[2]((1/4)*Pi) = 0}

Maple 2017 should be able to solve this set of linear differential equations as

a single whole system but it fails to do so when N > 2.  I don't know why.

Therefore I will guide it to arrive at a solution as follows.

The coupling of these differential equations is triangular in form; the solution may be

achieved by solving the equations in sequence, a single equation at a time.  This is

done in the following for-loop. The loop performs N assignments of the form:
    sol[i] := u[i](x) = an expression in x

where i=0,...,N-1.

for n from 1 to N do
  `if`(n=1, sys[1], eval(sys[n], {sol[i] $i=0..n-2}));
  dsolve(%);
  sol[n-1] := combine(%, trig);
end do:

Now assemble the pieces into the overall solution:

eval(S, {sol[i] $i=0..N-1});

z(x) = -(1/3)*cos(2*x)+(1/6-(8/45)*cos(x)*2^(1/2)-(1/90)*cos(4*x))*delta+(-(1/90)*cos(x)*2^(1/2)*Pi-(1/90)*cos(x)*2^(1/2)+(2/45)*sin(x)*2^(1/2)*x-(1/90)*2^(1/2)*cos(3*x)-(1/1050)*cos(6*x)+(7/270)*cos(2*x))*delta^2

 

 

Download perturbation2.mw

@Math Pi Euler My calculation shows that the leading term of the solution is −1/5 cos 2x.  The graph of tomleslie's  numerical solution, with its minimum at −0.2 (that is, −1/5) confirms that.

The leading term of your calculation, −1/3 cos 2x, disagrees with the evidence provided by the numerical solution.

Perhaps there is a typo in the differential equation that you gave us in your original post?

 

Which is the unknown that you want to solve for?

For everyone's convenience, here I transcribe Mathematica's solution to Burgers' equation into the Maple notation and verify that indeed it is correct.

restart;

A := 1 + erf(x/(2*sqrt(mu*t))):
B := 1 + erf((t-x)/(2*sqrt(mu*t))):
C := exp(-(t-2*x)/(4*mu)):
sol := u(x,t) = 1/(1 + C*A/B);

u(x, t) = 1/(1+exp(-(1/4)*(t-2*x)/mu)*(1+erf((1/2)*x/(mu*t)^(1/2)))/(1+erf((1/2)*(t-x)/(mu*t)^(1/2))))

pde := diff(u(x,t),t) + u(x,t)*diff(u(x,t),x) = mu*diff(u(x,t),x$2);

diff(u(x, t), t)+u(x, t)*(diff(u(x, t), x)) = mu*(diff(diff(u(x, t), x), x))

Verify that u(x, t) satisfies the PDE:

pdetest(sol, pde);

0

Verify that the initial condition on x > 0 holds:

limit(sol, t=0, right) assuming x>0, mu>0;

u(x, 0) = 0

Verify that the initial condition on x < 0 holds:

limit(sol, t=0, right) assuming x<0, mu>0;

u(x, 0) = 1

Download mw.mw

@Mariusz Iwaniuk Here Maple is really confused because u(x,t)=1 does not satisfy the initial condition.  As to the HINT=f(x)/f(t), that's a good try, but in view of the exact solution, we should be looking for HINT=f(x/t).  Unfortunately that does not result in any solution.

 

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);

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