Rouben Rostamian

MaplePrimes Activity


These are questions asked by Rouben Rostamian

A colleage asked me for the solution of an elementary boundary value problem for an ODE that models the steady-state temperature distribution in a nonhomogeneous rod.  I solved it in Maple and passed the solution on to her without checking the result.  She struggled for a day or two, attempting to "debug" her finite element code whose result was not agreeing with Maple's.

Upon closer inspection, it turned out that her code was correct and the solution returned by Maple was not.  Here are the details.  I can't tell whether Maple's dsolve can be improved to provide correct solutions to this and similar problems.

Problem with dsolve()

restart;

kernelopts(version);

`Maple 2025.0, X86 64 LINUX, Mar 24 2025, Build ID 1909157`

A boundary value problem for an ODE

 

Let's solve the following boundary value problem:

de := Diff(a(x)*Diff(u(x),x),x) = -1;

Diff(a(x)*(Diff(u(x), x)), x) = -1

bc := u(-1)=0, u(1)=0;

u(-1) = 0, u(1) = 0

Maple's dsolve() fails to find a solution:

dsolve({de,bc});

 

Solution obtained by hand

 

The solution may be obtained by hand, as follows.

 

Integrate the DE once

a(x)*(diff(u(x), x)) = c-x.

where c is the integration constant.  Therefore

diff(u(x), x) = (c-x)/a(x) and (c-x)/a(x) = c/a(x)-x/a(x).

To determine the value of the constant c, integrate the above over [-1, 1]

and apply the boundary conditions:

c*(int(1/a(x), x = -1 .. 1))+int(x/a(x), x = -1 .. 1) = u(1)-u(-1) and u(1)-u(-1) = 0

whence
c = (int(x/a(x), x = -1 .. 1))/(int(1/a(x), x = -1 .. 1)).

Having thus determined c, we integrate the expression for diff(u(x), x) obtained above, and arrive at the solution
u(x) = int((c-xi)/a(xi), xi = -1 .. x).

Remark: The solution obtained here is valid for any a(x) as long as

the integrals encountered above make sense.  Note that there is
no requirement on the differentiability or even continuity of a(x).

Technically, a(x) can be any function such that 1/a(x) is integrable.

 

 

Wrong solution returned by dsolve()

 

Let us consider a special choice of the coefficient a(x):

a := x -> 1 + 2*Heaviside(x);

proc (x) options operator, arrow; 1+2*Heaviside(x) end proc

According the the calculations above, we have:

c := int(x/a(x),x=-1..1) / int(1/a(x),x=-1..1) ;

-1/4

and therefore the correct solution is

int((c-xi)/a(xi), xi=-1..x):
sol := u(x) = collect(%, Heaviside);

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

Here is what the solution looks like:

plot(rhs(sol), x=-1..1, color=blue, thickness=3);

Let's attempt to solve the same boundary value problem with Maple's dsolve().

dsolve({de,bc}):
bad_sol := collect(%, Heaviside);

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

That's very different from the correct solution obtained above.   To see the problem with it, let's recall that according to the DE, the expression a(x)*(diff(u(x), x))NULLshould evaluate to  c-x for some constant c,  but what we get through bad_sol is nothing like c-x;  in fact, it's not even continuous:

a(x)*diff(rhs(bad_sol), x):
collect(%, Heaviside, simplify);
plot(%, x = -0.5 .. 0.5);

((4*x-2)*exp(-2)+4*x+2)*Heaviside(x)^2/(3*exp(-2)+3)+(-(4/3)*x-1/3)*Heaviside(x)-x-1/(3*exp(-2)+3)

It appears that dsolve goes astray by attempting to expand the expression diff(a(x)*(diff(u(x), x)), x).   It shouldn't.

 

 

Download dsolve-bug.mw

 

I seem to recall that a[k] and a(k) are treated interchangeably within a proc.  But I did not expect that to happen at the top level.  See the worksheet below.  Is there a way to keep indexed quantities as indexed?

restart;

kernelopts(version);

`Maple 2025.0, X86 64 LINUX, Mar 24 2025, Build ID 1909157`

A := piecewise(x<0, a[1], a[2]);

A := piecewise(x < 0, a[1], a[2])

int(1/A, x=-1..1);

1/a(1)+1/a(2)

Download mw.mw

I wish to plot two line segments over the intervals [0,1] and [2,3].  Maple 2024 draws a single line segment over [0,3].  Oddly enough, if I specify a color option, then it produces the expected result!

This used to work correctly in Maple 2023.  Perhaps this can be corrected in the forthcoming release?

restart;

kernelopts(version);

`Maple 2024.2, X86 64 LINUX, Oct 29 2024, Build ID 1872373`

plots:-display(
        plot(1, x=0..1),
         plot(1, x=2..3), thickness=5);

plots:-display(
        plot(1, x=0..1),
         plot(1, x=2..3),thickness=5, color=red);


Download mw.mw

restart;

Here are the graphs of a parabola and a straight line:

plots:-display(
        plot(x^2, x=-1..1),
        plot((x+1)/2, x=-1..1),
color=["Red","Green"]);

 

Suppose I want to plot the part of the parabola that lies below

the straight line, and suppose, just to be nasty, I choose to do it

with implicitplot:

plots:-implicitplot(y=x^2, x=-1..1, y=0..(x+1)/2);

 

That is not a parabola at all.  [And where does the "ynew" label come from?]

 

This behavior was introduced in Maple 2022.

In Maple 2021 we get the expected result:

plots:-implicitplot(y=x^2, x=-1..1, y=0..(x+1)/2);


 

Download mw.mw

 

restart;

Let x be some name:

x := asdf;

asdf

I wish to make a new name, y, whose value is the first character in x:

convert(x, string):
y := convert(%[1], name);

a

That works but seems too convoluted.  Is there a better way of doing that?

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