Rouben Rostamian

MaplePrimes Activity


These are replies submitted by Rouben Rostamian

@Markiyan Hirnyk Let me point out that the paragraph that yourself quoted from Wiki refers to the Dirac delta function.  In my post I used function in that same customary sense.  I trust that people who are knowledgeable about the subject are quite comfortable with that usage.

 

@Jjjones98 Here is how one derives the expression for the Green's function by hand.

Consider the fourth order ODE y''''(x) = Dirac(x-u). The Dirac function is zero to the left and to the right of the point x=u. Therefore we may solve the ODE y''''(x)=0 on the interval 0 < x < u and on u < x < 1 separately. Let's call those two solutions yL(x) and yR(x).

On the first interval we have:
yL''''(x)=0, yL(0)=0, yL'(0)=0, yL(u)=a, yL'(u)=b.

On the second interval we have:
yR''''(x)=0, yR(u)=c, yR'(u)=d, yR(1)=0, yR'(1)=0.

Here a, b, c, d are to be determined. Toward that end, you calculate yL(x) and yR(x) by hand. Then you determine the four constants a, b, c, d by solving the following system of four linear equations in the four unknowns a,b c, d:
  yL(u)=yR(u),
  yL'(u)=yR'(u),
  yL''(u)=yR''(u),
  yL'''(u) - yR'''(u) = 1.

Doing this is quite straightforward by hand, if you have the patience for it. Alternatively, just do it in Maple with dsolve() as I have shown. It does exactly the same thing.

@das1404 In my worksheet I refer to some equations through their labels.  That's not available in older Maples.  Here I am posting an edited version which avoids such references.  I saved the result as a *.mws file from within Maple 2017.  See if that works for you.

As to what you refer to as "Tom", I can't tell what it is intended to be since I am unable to view the result of your work.  When I execute your worksheet in Maple 2017, it stops with an error message before it hits the animation part.

Perhaps you can edit my worksheet and add the "Tom" part yourself.

Download cone-map.mws

@Jjjones98 One does not "prove" a definition.  One just "gives" a definition.

For instance, in your problem you should begin by defining what you mean be the Green's function.  Then investigate what that definition implies.  If you do things right, you will arrive the piecewise form that you gave in your original post.

You have said that you have proved various continuity and differentiability properties of your problem's Green's function.  What did you take as the definition of the Green's function?

Actually all the information you need is in Markiyan's post.  He gives a link to a wikipedia article where the Green's function is defined.  Then he uses Maple to derive the consequences of that definition.  His final result is not exactly what you are looking for but it takes just one small extra step to get there as shown here:

restart;

de := diff(y(x),x$4) = Dirac(x-u);

diff(diff(diff(diff(y(x), x), x), x), x) = Dirac(x-u)

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

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

sol := dsolve({de, bc}) assuming u > 0, u < 1;

y(x) = (1/2)*(u^3-2*u^2+u)*x^2-(1/2)*u*Heaviside(x-u)*x^2-(1/6)*Heaviside(x-u)*u^3+(1/2)*Heaviside(x-u)*u^2*x+(1/6)*(-2*u^3+3*u^2-1)*x^3+(1/6)*Heaviside(x-u)*x^3

simplify(convert(sol, piecewise, x));

y(x) = piecewise(x <= u, -(1/6)*x^2*(u-1)^2*(2*u*x-3*u+x), u < x, -(1/6)*u^2*(x-1)^2*(2*u*x+u-3*x))

 

Download mw.mw

@tomleslie That's very nice!

@Preben Alsholm Thanks for spotting and pointing out the typo in my worksheet.  I could see that the result was incorrect but was unable to locate the source of the error.

Additionally, thanks to torabi, who pointed out the typo in the definition of the function G.

After fixing those errors, the worksheet produces the figure shown in the cited article.

Here is the figure:

and here is the corrected worksheet: frac-dsolve-ver2.mw

 

@Daniel Skoog That's very nice.  I didn't know about that application.

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

 

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