mmcdara

7891 Reputation

22 Badges

9 years, 52 days

MaplePrimes Activity


These are replies submitted by mmcdara

@salim-barzani 

Remove the "&" character in the name of your file for it can be downloaded

and solve(%, R) returns 0 with multiplicity 2.
Explanation pde_mmcdara.mw

To complete @gkokovidis answer: as 0 is neither negative nor positve, it's worth looking to what happens when you used the Heaviside function:

restart
f := x -> Heaviside(-x)
x -> Heaviside(-x)
f~([-1, 0, 1])
                       [1, undefined, 0]
NumericEventHandler(invalid_operation = `Heaviside/EventHandler`(value_at_zero = 1)):
f~([-1, 0, 1])
                           [1, 1, 0]
NumericEventHandler(invalid_operation = `Heaviside/EventHandler`(value_at_zero = 1/2)):
f~([-1, 0, 1])
                           [   1   ]
                           [1, -, 0]
                           [   2   ]
NumericEventHandler(invalid_operation = `Heaviside/EventHandler`(value_at_zero = 0)):
f~([-1, 0, 1])
                           [1, 0, 0]
_EnvUseHeavisideAsUnitStep := true: 
f~([-1, 0, 1])
                           [1, 1, 0]

# Watch out: running _EnvUseHeavisideAsUnitStep := true:  
# makes `Heaviside/EventHandler`(value_at_zero = 0) ineffective
NumericEventHandler(invalid_operation = `Heaviside/EventHandler`(value_at_zero = 0)):
f~([-1, 0, 1])
                           [1, 1, 0]


Here are two tricks I use to ease the coding and avoid forgetting some backslashes (I used the substitute "&" in the second trick but you use any other character).
Based on @acer's answer:Sugggestions.mw

@acer 

I was completely off the mark in thinking that the solid line came from the way the procedure was called.
Now that you mention it, it seems to me that I've already noticed this effect due to the thickness of the line.

Thanks again for your comment and explanation.

@acer Thanks a lot.

I have two questions:

  1. As no Warning, `i` is implicitly declared local ...l is issued while i is not assigned to something (Maple 2015), why did you declare i as local? As a rule I explicitely declare as local all the names for which such a warning is fired: does that mean that the procedures I write this way are not safe?
     
  2. Could you explain me what happens in the las case
    ecdf(S, 'color'="green", 'thickness'=5)

    Why are the extreme segments in plain linestyle and not in dashed linestyle? I see that is related to the way you write the optional parameters but I still don't understand what happens here.

@dharr 

I have done a few attempts with "my" solution (assuming it  answer the OP's need, which I'm not sure):
 

initial_conditions := phi(0) = phi__0, D(phi)(0) = varphi__0;
sys  := {ODE, initial_conditions}:

Digits := 10:
sol := dsolve(
         sys
         , numeric
         , method=rkf45
         , events=[[phi(x)=0, halt]]
         , parameters=[phi__0, varphi__0]
       );


sol(parameters=[ two_numeric_values]);

interface(warnlevel=0): sol(1000):

XMAX := 1000:
xmax := sol(eventfired=[1])[];
if xmax=0. then xmax := XMAX end if:


DocumentTools:-Tabulate(
  [
    plots:-odeplot(sol, [x, phi(x)], x=0..xmax)
    ,
    plots:-odeplot(sol, [phi(x), v(phi(x))], x=0..xmax)
  ]
  , width=60
)


I was never capable to get plots closed to the the picture the OP gives.


Another test is based on a visual approximation of phi(x) from curve A:

test := 0.022*(1-tanh((x-150)/40))/2: 
plot(test, x=0..250);             # looks like curve A "phi(x) vs x"
plot([test, v(test), x=0..250])   # far from curve A "v(phi(x)) vs phi(x)"


Finally : v(Phi) is a decreasing negative valued function of Phi from Phi varying from 0 to 0.022 (ranges from curve A):

v(Phi);
plot(v(Phi), Phi=0..0.022)

So I can't figure out where the" Curve A pattern v(phi(x)) vs phi(x)" could come from?
 

@mary120 

Could you provide a file containing the values of x and phi(x) for curve A ?

Thanks

@salim-barzani 

In your initial qiuestion I simply followed step-by-step the paper you attached and it was clear and self consistent.
It's no longer the case here and I don't see the connection between this new paper and your attached worksheets.

For instance (file a[12]-pde.mw): you have now a nonlinear pde with integral terms whose unknown depends on 3 independent variables (a linear pde in x and t in the initial question).

  • How can you get rid of the nonlinear terms?
  • Do you want me to apply what I did in my answer to the linear part alone of this new equation?
  • What are the coefficients you want to equate to 0?
  • ... and so on.
     

Here is what you must take as my final contribution to this thread xxx.mw.
The procedure F generalizes procedure f in my initial answer, procedure `w'/w` has been update to handle more general situations.
I hope you will get more help for someone else.


By the way: Do you realize that you don't ask questions about Maple but that you ask helpers to do what is done in a paper, which requires a huge lot of time and enough knowledge of the underlying Physics? For you all this Physics is likely evidence but it's not the case for other Mapleprimes members, including me at first place.

 

@dharr 

... here is another point:

  • You surely noticed this "strange" IC
    initial_conditions := phi(0) = 0.022, D(phi)(0) = 0.0;

    and probably wondered why the OP wrote this curious D(phi)(0) = 0 extra condition?
    (Which is why you use in your worksheet 

    sol := dsolve([ODE2, initial_conditions[1]], phi(x), numeric, range = 0 .. 10, method = rkf45, output=listprocedure);

    )

  • I guess the OP doesn't try to solve ODE but diff(ODE, x): look to the patterns this could give:
     

    restart;

    # Define constants
    mu := 0.01:
    nu := 1 - mu:
    beta := 0.05:
    alpha := 0.3:
    M := sqrt(0.704):
    gamma1 := 0.001:

    # Define the function v(phi)

    v := f -> (1-alpha)*M^2 - (1-alpha)*M*sqrt(M^2 - 2*f) + mu*(mu + beta*nu)*(1 - exp(f/(mu + beta*nu))) +
                (nu/beta)*(mu + beta*nu)*(1 - exp(beta*f/(mu + beta*nu))) + (alpha/gamma1)*(1 - exp(-gamma1*f)):

    # Define the ODE (dphi/dx)^2 + v(phi) = 0

    ODE := diff(phi(x), x)^2 + v(phi(x)) = 0:
    ODE := diff(ODE, x):

    # Initial conditions: dphi/dx = 0 at x = 0, phi = 0 at x = 0
    initial_conditions := phi(0) = 0.022, D(phi)(0)=0.0007;

     

    phi(0) = 0.22e-1, (D(phi))(0) = 0.7e-3

    (1)

    sys  := {ODE, initial_conditions}:
    xmax := 100:
    sol  := dsolve(sys, numeric, range = 0..xmax, method = rkf45):

    plots:-odeplot(sol, [phi(x), v(phi(x))], x=0..xmax);

     

    plots:-odeplot(sol, [x, phi(x)], x=0..xmax);

     

     

    ``

  •  

    Download rk4_mmcdara.mw

     


Finally; are we sure the OP doesn't want to solve 

ODE := diff(phi(x), x$2) + v(phi(x));

?
In its initial problem v(phi(0)) is Imaginary, and then phi(x) is a complex valued function (whose numeric value along a given branch can be ontained using the dsolve option complex=true).
At no time the OP suggests the phi(x) is imaginary (it even writes "phi and x are real"), which makes me doubt about the correctness of the term diff(phi(x), x)^2

@salim-barzani 

It is pretty obvious that eqt3 does not verify pde.
Look the attached file for a numeric proof.
sol_mmcdara.mw 

# eqt verifies pde

simplify(eval(pde, eqt), size);
                               0

@dharr 

I wrote ":" instead of "-"

I edited my initial comment to show that choosing an a priori branch is not the good method.
In this edited worksheet no numeric soluion agrees with the one or the other of the two exact branches.
 

@acer 

I've carefully read your test cases number 9 to 11 and you raise questions I hadn't thought of, so I don't have any immediate answers to give you.
In fact (Test 9), I hadn't realized that there could be a difference between an “abstract” function and a predefined function in Maple (let's say that f(X) or cos(X) are, for me, the same thing).
I'm off to bed now and will get back to you tomorrow when I'm thinking more clearly.
Thanks again


@janhardo 
Thanks for the contribution.
I will look to your reply tomorrow, have a good day (night?).

@acer 

Thank you very much.

I verified the 7 test cases and, if I'm really fussy,only 3 of them H do not give me the expected answer.
It's probably about my forgetting that `*` and `^` are different and that I missed the 

  • or a linear combination of powers of (products of)  X1, ...XN . or functions of  X1, ...XN .

It could better if H( X1^3*X2^2,{X1, X2}) returned true, but as H(f(X1, X2), {X1, X1)}) = true, it is enough to define f that way
f := (u, v) -> u^3 * v^2 ... and everything's back to normal.

A priori, your procedure seems to meet all my expectations (provided possibly that the expression to be tested is written appropriately...but that's not a blocking point).

Here is my feedback: typefunc_ex_feedback.mw

PS I'm going to try and understand how your procedure works... but I've always thought that everything to do with properties and (structured) types is the most complex part of Maple and it's likely I will have some questions about it.

5 6 7 8 9 10 11 Last Page 7 of 154