mmcdara

7289 Reputation

18 Badges

8 years, 259 days

MaplePrimes Activity


These are replies submitted by mmcdara

@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?
Unless the initial differential equation bifurcates when some critical point is reached...

@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.

@janhardo 

I'd never thought of defining a new type this way. (I always use TypeTools[AddType]).
Thanks for the tip.

@nm 

Thanks for having shared your experience

@nm 

Thanks.

What do you usually use when you want to check if an expression has a given compound/structured type?

@nm 

I hadn't  thought of using the debugger.
Indeed an excellent idea, I'm going to try and see what I can get out of it.

Thank you very much.

@salim-barzani 

The "Big Sum" over i, j, ...m, n should be written as a summation over all n tuples (i, ..., n) which verify several conditions.
For whatever reason the authors of this paper have prefered not to account for these conditions to write this "Big Sum".
But the price to pay is that they have to divide its result by 2M*M! : the attached file (try to) explains where this corrective term comes from.

explanation.mw

@dharr 

Sure, you're right.
Thanks for the corrrection (and your vote)

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