Rouben Rostamian

MaplePrimes Activity


These are replies submitted by Rouben Rostamian

The article that you have cited is available for download from ResearchGate.

I spent some time trying to understand it, but was slowed down by the

sloppy presentation, inaccuracies, and omissions.  I had to figure out

what the author intends to say by examining his examples.  I must

say that I am not convinced of his method's validity.

 

In the first example, he solves the initial-boundary-value problem
for the unknown u(x, t)
"(&PartialD; u)/(&PartialD; t)= 1/(2)*x^(2)* ((&PartialD;)^2u)/((&PartialD;)^( )x^2) ",      " 0 < x < 1,     t >0,"

subject to the boundary conditions

"u(0,t)=0,     u(1,t)=e^(t),"
and the initial condition
u(x, 0) = x^2.

His method leads to a solution in the form of infinite series
"u(x,t)=(&sum;)`u__n`(x,t)=x^(2)*(1+t+1/(2!)*t^(2)+1/(3!)*t^(3)+...) = x^(2)*e^(t)"
which is the exact solution of the problem.  That's impressive.

 

But then, as an exercise, I applied his steps to solve the simpler problem
"(&PartialD; u)/(&PartialD; t)= ((&PartialD;)^2u)/((&PartialD;)^( )x^2) ,           0<x<1,    t>0, u(0,t)=0,    u(1,t)=0, u(x,0)=x*(1-x)."

In this case his series consists of a single term (the remaining terms

being all zeros)  leaving us with the "solution" u(x, t) = x*(1-x)

which is far from correct.

 

I may have misunderstood his method.  Try solving this problem with

his method and show me if you get something different.

 

An alternative to dharr's suggestion is to press Ctrl+Shift+G (that is, the Ctrl, Shift, and G keys together) followed by the Greek letter code.  For instance, Ctrl+Shift+G followed by e produces ε, and Ctrl+Shift+G followed by s produces σ.  Most Greek letters map to their corresponding Latin keys.  For the complete table of Greek character, look up the help page 2-D Math Shortcut Keys and Hints.

The Ctrl+Shift+G combination works on Linux and Windows.  On a Mac, type Command+Shift+G instead.

@Preben Alsholm What you have quoted can be an explanation of what we are seeing.  It's not quite clear to me how, but it seems relevant.  Thanks for your input.  Because of that I have converted your Reply to an Answer.

 

@Pascal4QM You need to say
A,C := test()[[1,3]];  

@mmcdara  Oh, I see what you mean—you want to consider the case where the beam is not mechanically attached to a support, and therefore it may detach and lift off.

That requires handling boundary condition given as inequalities rather than equalities.  That's a difficult problem, both mathematically and computationally.

However, there is some hope.  Beamsolve reduces Euler's PDE to a system of ODEs.  It is conceivable that the one-sided constraints at the supports may be handled through dsolve's Events mechanism as described in the help pages:

?dsolve,events
?How Do I Solve an Ordinary Differential Equation?

That can be a fun project for someone to take up.
 

@mmcdara Any of the values of u, uxuxx, and uxxx may be specified for the boundary conditions.  I have picked u=0 and uxx=0 as defaults because that corresponds to a simple support (a ball-and-socket support, as you put it) which I assume would be the most common type of support. I don't understand what you mean by "this prevents the simulation of simple supports", since the user's specification overrides the defaults.  Let me know if I have misunderstood you.

 

@rookie Beamsolve computes the values of u(x,t) in order to plot the solution but it does not provide a user-accessible function to extract those values.  I will think about adding one but this may take some time.  In fact, it will be good to have functions that return not only u, but also ux, uxx, and uxxx.

@mmcdara Thanks.  Calculating the stress is easy -- it's just a multiple of uxx.  Producing good-looking graphics to highlight it, is another story.  That can wait until after I am retired :-)

@Pemudahijrah01 Okay, then try this:
DEtools:-DEplot3d(sys, [x(t),y1(t),y2(t)], t=0..500, stepsize=0.5,
    [[x(0)=0.6, y1(0)=0.3, y2(0)=0.3]], arrows=cheap, linecolor=black);

@Pemudahijrah01 Is this good?
DEtools:-DEplot(sys, [x(t),y1(t),y2(t)], t=0..50, stepsize=0.1,
    [[x(0)=0.6, y1(0)=0.3, y2(0)=0.3]], scene=[t,x(t)], x=0..0.6,
    linecolor=black, thickness=1);

@Pemudahijrah01 When you change Mu[2]=0.74 to Mu[2]=0.38, the equilibrium shifts to  {x = 0.5233951682, y1 = 0.2697791086, y2 = 0.}.  Note the y2=0 part!

The new equilibrium is still stable.  The x(t) and y1(t) converge quickly to their equilibrium values, but y2(t) takes a very long time to go there.

To adjust the phase portrait for the new parameter, change the previous
y2min, y2max := 0.1, 0.2;
to
y2min, y2max := 0.0, 0.01;
Then we get:

It helps to plot x(t), y1(t), y2(t) as functions of time.  These show that x(t) and y1(t) stabilize within 0 < t < 50, but for y2(t) you need 0 < t < 1000.

DEtools:-DEplot(sys, [x(t),y1(t),y2(t)], t=0..50, stepsize=0.1,
    [[x(0)=xmax, y1(0)=y1max, y2(0)=y2max/2]],
    x=xmin..xmax, y1=y1min..y1max, y2=y2min..y2max, scene=[t,x(t)],
    linecolor=black, thickness=1);

@rookie Here is a version of the worksheet in which I have replaced some of the newer Maple syntax with their older equivalents.  This one works with Maple 2017 (tested) and perhaps earlier (not tested).

Download worksheet: euler-beam-with-method-of-lines-2017.mw

 

Here are a couple of things to think about:

  1. Your data set consists of 42 points.  Simpon's rule of integration requires an odd number of data points. Why are you focused on Simpson's rule in this case?
  2. The factor sin(x) looks odd.  The "angle" x goes from 1 to 42.  What are the units of measurement here? Does sin(42) mean the sine of 42 degrees, 42 grads, or 42 radians?  Unless you qualify this, Maple will assume that you mean 42 radians.  Is that what you intend?

 

@mary120 The "trick", if it can be called so, is to plot the system's vector field first, without worrying about the solution curves.  Thus, you will do:
DEtools:-DEplot({de1,de2}, [x(t),y(t)], t=0..1, x=-4..4, y=-4..4);
and get:

From that picture you can see that if you start anywhere on the positive x axis, let's say [x(0)=3, y(0)=0], as you have done, the orbit (that's the correct word for what you have called "solution curves") takes off from that point and moves up and to the right until it exits the right edge of the plot.  Something that may not be immediately obvious to you is that if you start at that same initial point and move backward in time, that is, go against the arrows, you will move down and to the right until you exit the right edge of the plot.  Be sure to understand what I am saying here.

Therefore, to sketch the orbit starting at [x(0)=3, y(0)=0], going both backward and forward in time (note the t=−1..1), we do

DEtools:-DEplot({de1,de2}, [x(t),y(t)], t=-1..1, [[x(0)=3, y(0)=0]], x=-4..4, y=-4..4, stepsize=0.01);

Once you see this, you can add a sequence of initial conditions spread along the x axis through the seq() constructions which I showed in my previous answer.  To complete the diagram, you will then add initial conditions along the y axis in the same way.  Then you will adjust the line thickness and color, again as I did in my previous answer.

If you want to ask further questions, include your worksheet to show what you have done.

 

 

 

@dharr I agree with your observation regarding user-friendliness of error messages.  As much as it is tempting to perform type-checking at the level of parsing a proc's args, that can lead to cryptic error messages as you have noted.

Customized parsing of the args adds a few extra lines of code but that can be worth the effort in a proc which is meant to be called at the top-level by other users.  And as you have noted, that "other user" can be you!

That said, Carl in his responses has pointed out that there are cases where type-checking at the argument level can be justified.  That's something to have in mind when writing a general-purpose code that's expected to have a lifetime of more that a few days.

 

 

First 37 38 39 40 41 42 43 Last Page 39 of 99