Rouben Rostamian

MaplePrimes Activity


These are replies submitted by Rouben Rostamian

@Kitonum That's nice and works well with low degree polynomials.  In the general case, I would feel uncomfortable in not knowing how Maple is going to parametrize the implicitly defined function, or if it is going to be able to parametrize it at all.  The direct differentiation method would be a better approach in those cases.

@Jesús Guillera and @Kitonum, here is the correct calculation of the second derivative.  Sorry for the previous misfire.

restart;

P:=(x,y)->11*x^3-9*x^2*y+27*x*y^2+7*y^3-18*x^2-24*x*y+18*y^2;

proc (x, y) options operator, arrow; 11*x^3-9*x^2*y+27*x*y^2+7*y^3-18*x^2-24*y*x+18*y^2 end proc

plots:-implicitplot(P(x,y), x=-1..3, y=-2..3, gridrefine=3, scaling=constrained, color=black, thickness=3);

The first derivative

 

Consider y as a function of x and calculate the derivative

P(x,y(x));
diff(%, x):
dP := simplify(isolate(%, diff(y(x), x)));

11*x^3-9*x^2*y(x)+27*x*y(x)^2+7*y(x)^3-18*x^2-24*y(x)*x+18*y(x)^2

diff(y(x), x) = (9*y(x)^2+(-6*x-8)*y(x)+11*x^2-12*x)/(-7*y(x)^2+(-18*x-12)*y(x)+3*x^2+8*x)

Introduce notation for the numerator and denominator of that derivative:

nn := numer(rhs(dP));
dd := denom(rhs(dP));

-11*x^2+6*y(x)*x-9*y(x)^2+12*x+8*y(x)

-3*x^2+18*y(x)*x+7*y(x)^2-8*x+12*y(x)

We want to evaluate the derivative at the origin, but the numerator and denominator

are zero there:

subs(x=0, y(0)=0, nn);
subs(x=0, y(0)=0, dd);

0

0

Apply l'Hopital's rule to resolve the indeterminacy:

diff(nn, x) / diff(dd, x);
limit(%, x=0);
L := simplify(subs(y(0)=0, %));

(-22*x+6*(diff(y(x), x))*x+6*y(x)-18*y(x)*(diff(y(x), x))+12+8*(diff(y(x), x)))/(-6*x+18*(diff(y(x), x))*x+18*y(x)+14*y(x)*(diff(y(x), x))-8+12*(diff(y(x), x)))

-(9*y(0)*(D(y))(0)-3*y(0)-4*(D(y))(0)-6)/(7*y(0)*(D(y))(0)+9*y(0)+6*(D(y))(0)-4)

(3+2*(D(y))(0))/(-2+3*(D(y))(0))

That should be the same as the derivative D(y)(0).  Therefore we equate and solve:

L = D(y)(0):
S := solve(%, D(y)(0));

2/3-(1/3)*13^(1/2), 2/3+(1/3)*13^(1/2)

Those are the slopes of the two branches of the curve at the origin.
Their numerical values agree with those calculated by Jesús Guillera:

evalf([S]);

[-.5351837583, 1.868517092]

 

The second derivative

 

We are going to calculate the second derivatives of the two branches,

based on the following data:

data1 := x=0, y(0)=0, D(y)(0)=S[1];
data2 := x=0, y(0)=0, D(y)(0)=S[2];

x = 0, y(0) = 0, (D(y))(0) = 2/3-(1/3)*13^(1/2)

x = 0, y(0) = 0, (D(y))(0) = 2/3+(1/3)*13^(1/2)

Begin with calculating the general form of the second derivative

diff(P(x,y(x)), x, x):
simplify(solve(%, diff(y(x), x, x))):
d2P := convert(%, D);

((18*x+14*y(x)+12)*(D(y))(x)^2+(-12*x+36*y(x)-16)*(D(y))(x)+22*x-6*y(x)-12)/(-7*y(x)^2+(-18*x-12)*y(x)+3*x^2+8*x)

As in the previous section, the numerator and denominator vanish at the origin

nn := numer(d2P):
dd := denom(d2P):

simplify(subs(data1, nn));
simplify(subs(data2, nn));
subs(x=0, y(0)=0, dd);

0

0

0

So we apply l'Hopital.  Here is the second derivative of one branch

limit(diff(nn, x) / diff(dd, x), x=0):
simplify(subs(data1, %)) = (D@@2)(y)(0):
solve(%, (D@@2)(y)(0));
evalf(%);

-209/81+(1057/1053)*13^(1/2)

1.039000660

and here is the second derivative of the other branch:

limit(diff(nn, x) / diff(dd, x), x=0):
simplify(subs(data2, %)) = (D@@2)(y)(0):
solve(%, (D@@2)(y)(0));
evalf(%);

-209/81-(1057/1053)*13^(1/2)

-6.199494488

 

 

Download second-derivative.mw

@Kitonum Yours is a good example and it shows that there is no basis to my previous calculation.  The calculation would have been correct if the expression returned by implicitdiff were continuous at (u0,v0) but it certainly is not, and therefore taking the limit toward (u0,v0) is not meaningful.

If I come up with a valid argument, I will write again.  In the meantime, thanks for your insightful comment.

You write: "We see that P(u0,v0)=0."

I don't see that.  To get more help, upload your worksheet by clicking the big fat green arrow when you reply to this message.

@janhardo As you have noted, your code needs several corrections.  But even after those corrections a major impediment remains.  Your code intends to apply pdsolve,numeric to find w(x,y,t).  The problem that we face is that Maple's pdsolve,numeric is designed to solve PDEs in two independent variables.  So we are out of luck with solving for w(x,y,t).

Yet another major issue: After a few corrections, the PDE that you have attempted to present would describe the equation of motion of a membrane, not a plate which is the subject of this post.  The equation of the plate is of 4th order in the spatial variables.  See the equation included in my earlier answer.

In your description of the boundary conditions, you refer to the "clamped" boundary condition.  There is no such a thing as a clamped boundary condition for a membrane.  The clamped boundary condition applies to plates, and indicates that not only the displacement, but also the slope at the boundary is zero. 

@Scot Gould Your PDE problem is quite well-posed.  It's not your fault that Maple is returning the wrong solutions.  That's a Maple bug.  See the stripped-down example that I posted a few minutes ago that comes with an exact solution, while Maple gives the wrong answer.

So, for now, you need to solve your prolem in the old-fashioned way by hand.  Hopefully Maple's pdsolve will be corrected in future releases.

@dharr The heat equation with any type of Robin boundary condition is always well-posed.  Physical considerations do not enter the picture.  I suspect that your remarks stem from the unstated desire to have solutions that remain bounded in time.  But that's not a requirement for the solvability of a PDE.

The following example shows that the boundary conditions of the sort that was originally posed by Scott Gould is quite viable.  The solution is easily constructed by hand.  Maple 2024.1 returns the wrong solution.

restart;

pde := diff(u(x,t),t) = diff(u(x,t),x,x);

diff(u(x, t), t) = diff(diff(u(x, t), x), x)

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

u(0, t)+(D[1](u))(0, t) = 0, u(1, t)+(D[1](u))(1, t) = 0

ic := u(x,0) = exp(-x);

u(x, 0) = exp(-x)

pdsolve({pde, bc, ic});

u(x, t) = 0

That is obviously wrong since it does not satisfy the initial condition.

 

But the solution can be computed by hand through a straighforward separation

of variables technique, whereby we obtain

ans := u(x,t) = exp(t-x);

u(x, t) = exp(t-x)

pdetest(ans, [pde, bc, ic]);

[0, 0, 0, 0]


Download robin.mw

@Gabriel Barcellos If your goal is to plot m versus T and G versus T, you don't need to calculate a sequence of points first.  You can plot them directly.  Here is how.

restart;

with(plots):

plots[setoptions](font=[TIMES,18], labelfont=[TIMES,24]);

eq1 := m = tanh(6*m/T);

m = tanh(6*m/T)

eq2 := G = -T*ln(2*cosh(6*m/T)) + 3*m^2;

G = -T*ln(2*cosh(6*m/T))+3*m^2

The graph of m versus T

 

Equation eq1 provides a relationship between m and T.

Certainly m=0 is a solution for any T>0.

But if T<6, we have two additional solutions of the form m=+a and m=-a

where a depends on T.  The following graph of m versus T illustrates that:

implicitplot(eq1, T=0..9, m=-1..1,
        color=red, thickness=5, view=[0..9, -1.2..1.2]);

Further down, we will need the equation of the curved part of

that graph.  The graph is symmetric about the horizontal axis,

so we need the equation of the upper half of that curve only.

The following proc receives an argument T < 6 and returns the
value of the corresponding m value.
As the curve does not go beyond T=6, the proc fails

to return a value if T>=6. That's by design.

F := proc(T)
        if not type (T, realcons) then return 'procname'(args) end if;
        fsolve(m = tanh(6*m/T), m=1e-9..1);
end proc:

For example, we have

F(4);

.8585596366

 

The graph of G versus T

   

Eq2 gives G as a function of m and T.  But m = F(T) as we saw

in the previous section. Therefore G is determined by T alone.

We wish to plot G versus T, but let us first
make a couple of observations.

 

1. m=0 is a solution of eq1 for all T.  Therefore we plug

m=0 in eq2 and obtain

    G = -T ln(2).

The graph of this branch of the solution is just a straight line.

 

2. If 0<T<6 we have two more solutions to eq1 given by m = F(T)

and m=-F(T).  But G is an even function of m, and therefore both

these solutions yield the same value of G.

 

We conclude that G versus T consists of two graphs.  One extends

for all values of T>0, and the other only for 0<T<6.  Here's what

those graphs look like.

eval(eq2, m=0);
p1 := plot(rhs(%), T=0..9, color="Red", thickness=5);

G = -T*ln(2)

myG := subs(m=F(T), eq2);
p2 := plot(rhs(%), T=0..6, color="Green", thickness=5);

G = -T*ln(2*cosh(6*F(T)/T))+3*F(T)^2

 

 

The graph of G versus T

display(p1, p2, labels=[T,G]);

Download mw.mw

 

 

I was bothered by the fact that the simple answer, b=11, was obtained after combining some considerably large and messy expressions.  I looked and found a very simple solution that avoids the mess.  Maple is not of much help in this alternative approach, so I will just post my non-Maple calculation here.

@vv That's much better than my brute-force calculation, and also more correct in accounting for the order of tangency.  Vote up!

Your equation is y' = (1-x)*y^(-1).  That's the Benoulli equation with P(x)=0, Q(x)=1-x, n=-1.

@imparter 

There have been over 30 exchanges in this thread, and we are not yet close to obtaining the results that you are looking for. It seems to me that you are attempting to program in Maple at a level which is beyond your current state of knowledge. Perhaps you should set this project aside and focus on learning the basics of Maple programming, and then tackle this more advanced problem after gaining sufficient mastery of the basics.  Attempting to force your way through solving this particular problem is not a good way of learning Maple.

@imparter A few more errors fixed.  Pay closer attention to the details.

Query-Comments2.mw

@imparter In the attached worksheet I have pointed out a few issues that require your attention.  See if you can fixed those.

Query-Comments.mw

@imparter This worksheet handles the two sets of data through a for-loop.

Two_data-Fixed.mw

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