Rouben Rostamian

MaplePrimes Activity


These are replies submitted by Rouben Rostamian

@vv Certainly the u(x,t)=sin(x)*exp(t) solves the backward heat equation exactly.  However it walks a very fine line.  Nearby functions generally evolve quite differently and that becomes a major issue under numerical noise.  Change t=0..1 to t=0..3 in your example to see that. We can hardly say that the result "fully agrees with the exact solution".

 

@Carl Love Thanks for pointing to the OP= help page.  I had not noticed this extremely helpful addition in Maple2019.

As to the wrong sign in the second derivative, yes, it is a mathematical flaw.
To understand that, we need some minimal information about the Sobolev
space 1/H(a, b).

The Sobolev space 1/H(a, b)contains just about any function and function-like
object that we can ordinarily imagine on the interval a, b, including differentiable
functions, discontinuous functions, and even things like the Dirac delta function
(which are not functions at all in the usual sense.)

It can be shown that the initial-boundary value problem for the heat equation
"(∂ u)/(∂ t)=((∂)^2u)/((∂)^( )x^2),                    x  in (a,b),  t>0,  u(a,t)=u(b,t)=0,    t>0,  u(x,0)=`u__0`(x),            "
`in`(x, a, b),
is solvable for any "`u__0` in "1/H(a, b).  The surprising thing is that the solution
at any t > 0  is in `#msup(mi("C"),mo("∞"))`(a, b).  Thus, the heat equation smoothes out any initial
condition in 1/Hto a C^infinity function instantaneously!.

Now, imagine that you reverse the direction of time, that is, you begin with a
function u__0 and you wish to evolve it backward in time via the heat equation.
What I have said above implies that  even if you begin with the nicest `in`(u__0, C^infinity),
the solution will lose its smoothness instantaneously and turns into things
like Dirac's delta or worse.  A numerical solver cannot deal with that.

Going backward in time via the heat equation amounts to changing t to -t,
whereby the equation turns into
"(∂ u)/(∂ t)=-((∂)^2u)/((∂)^( )x^2)".    (note the minus sign!)
The original PDE that was posed at the beginning of this thread was a more
complex version of this equation, that's why I commented that the sign of
the second derivative term is incorrect.

@arshl Solving your modified problem leads to an integral of the type
which is not expressible in terms of elementary functions. This means that your problem has no symbolic solution.  If this problem arises in a real application, then you may seek a numerical solution.  But for that you need to supply the numerical values of all the parameters.

 

 

@Carl Love Your approach requires the initial condition to be twice differentiable.  We can see that by replacing the initial condition u(x,0)=x*x(x-1) by u(x,0)=g(x) and then observing the resulting expression for bounds.

Thus, for instance, this approach fails with the perfectly fine initial condition u(x,0)=Heaviside(x-1/2).

If, however, the twice differentiability condition on the initial condition is met, then your approach will produce the correct solution.

That said, let me reiterate that the original statement of the PDE is flawed.  The sign of the second derivative term is wrong.  To see the problem, in your solution change
psol:= pdsolve(pde, bounds, numeric):
to
psol:= pdsolve(pde, bounds, numeric, spacestep=s);
and watch what happens to the amplitude of the solution with s=0.01 and s=0.001.

Question: The use of 'union=' in your solution is new to me. Where is it documented?

Added later:
I see that C-style assignments such as
x := 12;
x += 1;
x *= 5;

also work.  These have been on my wish-list for a long time!

@abbastalebi Here is how:

plots:-odeplot(dsol, [eta, diff(f(eta),eta)]);      # first derivative
plots:-odeplot(dsol, [eta, diff(f(eta),eta,eta)]);  # second derivative

@Mariusz Iwaniuk That's excellent. A thumb up!  I wonder whether it's possible to recover the original integral through differentiating the answer.

@tomleslie No, I didn't mean to.  I played around with the equation a bit and, among other things, changed the sign of the exponent, but then I forgot to restore it.  Mea culpa!

That said, the function exp(t*x) grows astronomically large very quickly.  I can't think of any realistic problem leading to such a term.  I like the problem better with the minus in the exponent.

@kambiz1199 A complete answer to your question depends on the details of your worksheet that you have not shown, so people are just guessing what you may have in there.

Here is a complete worksheet where for simplicity I have taken the N matrix to be 2x4, and have taken some random expressions for N[1] and N[2].  If this works for you, then you should be able to extend it as needed.

restart;

interface(rtablesize=20):

N[1] :=  x^2 + y^2;
N[2] := x^3*y + 7*y^3;

x^2+y^2

x^3*y+7*y^3

Nmat := < seq('N[i],0', i=1..2);  
          seq('0,N[i]', i=1..2) >;

Matrix(2, 4, {(1, 1) = x^2+y^2, (1, 2) = 0, (1, 3) = x^3*y+7*y^3, (1, 4) = 0, (2, 1) = 0, (2, 2) = x^2+y^2, (2, 3) = 0, (2, 4) = x^3*y+7*y^3})

Z := < D[1], 0; 0, D[2]; D[1], D[2] >;

Matrix(3, 2, {(1, 1) = D[1], (1, 2) = 0, (2, 1) = 0, (2, 2) = D[2], (3, 1) = D[1], (3, 2) = D[2]})

Matrix(3,4, (i,j)-> add(Z[i,k](unapply(Nmat[k,j], [x,y]))(x,y),k=1..2));

Matrix(3, 4, {(1, 1) = 2*x, (1, 2) = 0, (1, 3) = 3*x^2*y, (1, 4) = 0, (2, 1) = 0, (2, 2) = 2*y, (2, 3) = 0, (2, 4) = x^3+21*y^2, (3, 1) = 2*x, (3, 2) = 2*y, (3, 3) = 3*x^2*y, (3, 4) = x^3+21*y^2})


 

Download mw2.mw

@arshl  We apply dsolve to find the solution v[3](t) of the differential equation.  The solution depends on many parameters such as p1, R, k11, a11, ro, cp2, etc.  There is nothing special about p1.  Maple is quite happy solving differential equations involving parameters.  Here is a simple example:

de := diff(x(t),t,t) + omega^2 * x(t) = 0;
dsolve(de, x(t));

The solution x is a function of t, so we write it as x(t).  It also depends on the parameter omega, but we don't write x(t,omega).

 

@tomleslie Thanks for pointing this out.  I should have noticed that the integral symbol returned by dsolve is colored gray, indicating an inert Int.

 

@arshl In my previous reply I said that all "pi" should be "Pi".  You haven't made that change yet.  You should.

Looking further down your code, I don't understand the meaning of 
m1(p1, my(t)) := (1/2)*b11*ro*pi*R^2*(my(t)*cp2/(R*p1)-cp3)/p1;
What are your attempting to do there?

Perhaps it would be easier if you explain just the mathematics of your problem (no Maple).  Then someone may suggest ways to convert that mathematics into Maple.

Your say you are applying dsolve, but you haven't shown the code that produces the result that you have shown.  Showing only a part of your worksheet is not very informative.  Post your entire worksheet as an attachment.  How to do that? In the window in which you edit your message before posting, note the big green arrow near the top edge.  Click on that arrow to attach your worksheet.

But...

Before you do that, change all occurrences of "pi" in your worksheet to "Pi".  They are not the same thing.

 

@escorpsy

restart;

We have seen that the integral

exp(-h*Pi*v*(x+y)/beta) / cosh(Pi*v*(y-x)/beta)^(4+h):
Int(%, y=-infinity..z, x=-infinity..w);

Int(exp(-h*Pi*v*(x+y)/beta)/cosh(Pi*v*(y-x)/beta)^(4+h), y = -infinity .. z, x = -infinity .. w)

may be evaluated in terms of elementary functions when h is a given
negative integer.  In your latest worksheet you ask whether it is possible
to obtain the answer in terms of an unspecified h.

 

That is doubtful.  Consider the much simpler case

int(1/cosh(x)^h, x);

int(1/cosh(x)^h, x)

showing that Maple is unable to obtain an explicit formula for general h.

 

Specifying h as a non-integer does not help either:

int(1/cosh(x)^(7/3), x);

int(1/cosh(x)^(7/3), x)

But when h is an integer, we get an explicit formula:

int(1/cosh(x)^h, x) assuming h::posint;

(Sum((Product(1-1/(-h+2*j+1), j = 0 .. i))*sech(x)^(h-2*i-2)/h, i = 0 .. -ceil(-(1/2)*h)-1))*tanh(x)+2*(Product(1+1/(-h+2*j-1), j = 1 .. -ceil(-(1/2)*h)))*arctan(exp(x))

Unfortunately the result is too complex for evaluating the integral on -infinity, 0.

int(1/cosh(x)^h, x=-infinity..0) assuming h::posint;

limit(-(Sum((Product((h-2*j)/(h-2*j-1), j = 0 .. i))*sech(x)^(h-2*i-2)/h, i = 0 .. -ceil(-(1/2)*h)-1))*tanh(x)-2*(Product((h-2*j)/(h-2*j+1), j = 1 .. -ceil(-(1/2)*h)))*arctan(exp(x))+(1/2)*(Product((h-2*j)/(h-2*j+1), j = 1 .. -ceil(-(1/2)*h)))*Pi, x = -infinity)

Your actual integral is more complex than the simple illustration above,

therefore I wouldn't hold out hope for an explicit formula for its value.

 

Download mw.mw

@mmcdara I pointed out the ambiguity of the definitions of c_1 and c_2, but farazhedayati said that he has reasons for treating xbar as an independent variable, so I just implemented what he is asking.  It's up to him to make sense of the result.

 

@farazhedayati To add the c_2 values, we do

add(c_2(my_x, my_xbar, i), i=1..n);

and we get 1.479074043.

First 48 49 50 51 52 53 54 Last Page 50 of 99