mmcdara

6715 Reputation

18 Badges

8 years, 163 days

MaplePrimes Activity


These are replies submitted by mmcdara

@AHSAN 

restart

``

P0 := -6*k/((k+1)*(k-1)*(k*x-k-x)^2)-6/((k-1)*(k*x-k-x))-6/(k^2-1):

lambda[0] := k/(k+1):

P1 := -6*x*((x-1)^4*(A+(-1/3-4*Br*(1/21))*m)*k^8-(1/5)*(4*(((A*beta+15*Br*m*(1/56))*lambda[0]^3+(-9*A*beta*(1/5)-15*Br*m*(1/28))*lambda[0]^2+4*A*beta*lambda[0]*(1/3)+(-7*beta*(1/18)+5)*A-(1/21)*(20*(Br+7/4))*m)*x+(-2*A*beta-15*Br*m*(1/28))*lambda[0]^3+(18*A*beta*(1/5)+15*Br*m*(1/14))*lambda[0]^2-8*A*beta*lambda[0]*(1/3)+7*A*beta*(1/9)))*(x-1)^3*k^7+12*(x-1)^2*(((A*beta+15*Br*m*(1/56))*lambda[0]^3+(-9*A*beta*(1/5)-5*Br*m*(1/7))*lambda[0]^2+4*A*beta*lambda[0]*(1/3)+(-14*beta*(1/27)+5/2)*A-(1/21)*(10*(Br+7/4))*m)*x^2+((-A*beta-15*Br*m*(1/56))*lambda[0]^3+(9*A*beta*(1/5)+15*Br*m*(1/14))*lambda[0]^2-4*A*beta*lambda[0]*(1/3)+7*A*beta*(1/9))*x+(-4*A*beta*(1/3)-5*Br*m*(1/14))*lambda[0]^3+(12*A*beta*(1/5)+5*Br*m*(1/28))*lambda[0]^2-16*A*beta*lambda[0]*(1/9)+7*A*beta*(1/54))*k^6*(1/5)-(1/5)*(12*(((A*beta+5*Br*m*(1/28))*lambda[0]^3+(-9*A*beta*(1/5)-15*Br*m*(1/14))*lambda[0]^2+8*A*beta*lambda[0]*(1/9)+(-7*beta*(1/9)+5/3)*A-(1/63)*(20*(Br+7/4))*m)*x^3+((5/14)*Br*m*lambda[0]^3+(15/14)*Br*m*lambda[0]^2+(16/9)*A*beta*lambda[0]+(7/9)*A*beta)*x^2+((-A*beta-45*Br*m*(1/56))*lambda[0]^3+(9*A*beta*(1/5)+15*Br*m*(1/28))*lambda[0]^2-4*A*beta*lambda[0]+7*A*beta*(1/18))*x-2*lambda[0]*((A*beta+5*Br*m*(1/56))*lambda[0]^2-9*A*beta*lambda[0]*(1/5)+4*A*beta*(1/9))))*(x-1)*k^5+(((4*A*beta*(1/5)-3*Br*m*(1/7))*lambda[0]^3+(-72*A*beta*(1/25)-12*Br*m*(1/7))*lambda[0]^2-32*A*beta*lambda[0]*(1/15)+(1-56*beta*(1/45))*A-(1/21)*(4*(Br+7/4))*m)*x^4+((15*Br*m*(1/7)+4*A*beta*(1/5))*lambda[0]^3+(6*Br*m*(1/7)+144*A*beta*(1/25))*lambda[0]^2+32*A*beta*lambda[0]*(1/3)+28*A*beta*(1/45))*x^3+((-15*Br*m*(1/14)+4*A*beta*(1/5))*lambda[0]^3+(-396*A*beta*(1/25)+9*Br*m*(1/7))*lambda[0]^2-16*A*beta*lambda[0]*(1/3)+14*A*beta*(1/15))*x^2+(1/5)*(4*((-75*Br*m*(1/56)+A*beta)*lambda[0]^2+81*A*beta*lambda[0]*(1/5)-20*A*beta*(1/3)))*lambda[0]*x-32*A*lambda[0]^2*(-27/40+lambda[0])*beta*(1/5))*k^4+(((4*A*beta*(1/5)+9*Br*m*(1/14))*lambda[0]^3+(3*Br*m*(1/7)+108*A*beta*(1/25))*lambda[0]^2+16*A*beta*lambda[0]*(1/5)+14*A*beta*(1/45))*x^4+((-4*A*beta-9*Br*m*(1/14))*lambda[0]^3+(-324*A*beta*(1/25)+3*Br*m*(1/7))*lambda[0]^2-16*A*beta*lambda[0]*(1/5)+14*A*beta*(1/45))*x^3+8*lambda[0]*((A*beta-3*Br*m*(1/28))*lambda[0]^2+27*A*beta*lambda[0]*(1/25)-8*A*beta*(1/15))*x^2-8*A*lambda[0]^2*(lambda[0]-27/25)*beta*x-16*A*beta*lambda[0]^3*(1/5))*k^3-12*x*lambda[0]*(((A*beta+5*Br*m*(1/56))*lambda[0]^2+9*A*beta*lambda[0]*(1/5)+4*A*beta*(1/9))*x^3+((-3*A*beta+5*Br*m*(1/56))*lambda[0]^2-9*A*beta*lambda[0]*(1/5)+4*A*beta*(1/9))*x^2+2*A*lambda[0]*(lambda[0]-6/5)*beta*x+2*A*beta*lambda[0]^2)*k^2*(1/5)+12*x^2*A*((lambda[0]+3/5)*x^2+(-lambda[0]+3/5)*x-4*lambda[0]*(1/3))*lambda[0]^2*beta*k*(1/5)-4*A*x^3*beta*lambda[0]^3*(x+1)*(1/5))*(x-1)*(k-1)/(k^4*(k+1)*((x-1)*k-x)^6):

PP := P1*epsilon+P0:

ind := [indets(PP, name)[]];

[A, Br, beta, epsilon, k, m, x]

(1)

ind := convert(indets(PP) minus {x}, list);

epsilon_values := [.1, .3, .5, .7, .9];
k_values       := [seq(i, i = 2.5 .. 20, .5)];

data := (eval, kval) -> ind =~ [0.5, 0.2, 0.5, eval, kval, 1.5]:

# example

data(.2, 2.5)

[A, Br, beta, epsilon, k, m]

 

[.1, .3, .5, .7, .9]

 

[2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10.0, 10.5, 11.0, 11.5, 12.0, 12.5, 13.0, 13.5, 14.0, 14.5, 15.0, 15.5, 16.0, 16.5, 17.0, 17.5, 18.0, 18.5, 19.0, 19.5, 20.0]

 

[A = .5, Br = .2, beta = .5, epsilon = .2, k = 2.5, m = 1.5]

(2)

KE     := numelems(k_values)*numelems(epsilon_values);
AllMax := Matrix(KE, 3):

counter := 1:
for varepsilon in epsilon_values do
  for kappa in k_values do
    __PP             := eval(PP, data(varepsilon, kappa));
    Pre_prime        := diff(__PP, x);
    Pre_double_prime := diff(Pre_prime, x);
    critical_points  := solve(Pre_prime = 0, x, real, maxsols = 4);

    local_max := -infinity:
    for cp in [critical_points] do
      if eval(Pre_double_prime, x=cp) < 0 then
        local_max := max(local_max, eval(__PP, x=cp))
      end if:
    end do:
    AllMax[counter, 1] := varepsilon:
    AllMax[counter, 2] := kappa:
    AllMax[counter, 3] := local_max:
    counter := counter+1:
  end do:
end do;

AllMax;

Statistics:-ScatterPlot3D(AllMax, labels=[``(epsilon), ``(k), ``])

KE := 180

 

Matrix(%id = 18446744078491512094)

 

 

col := 1;
map(ke -> if AllMax[ke, col]=0.3 then AllMax[ke, [3-col, 3]] end if, [$1..KE]);

col := 1

 

[Vector[row](%id = 18446744078473882910), Vector[row](%id = 18446744078473883030), Vector[row](%id = 18446744078473883150), Vector[row](%id = 18446744078473883270), Vector[row](%id = 18446744078473883390), Vector[row](%id = 18446744078473883510), Vector[row](%id = 18446744078473871358), Vector[row](%id = 18446744078473871478), Vector[row](%id = 18446744078473871598), Vector[row](%id = 18446744078473871718), Vector[row](%id = 18446744078473871838), Vector[row](%id = 18446744078473871958), Vector[row](%id = 18446744078473872078), Vector[row](%id = 18446744078473872198), Vector[row](%id = 18446744078473872318), Vector[row](%id = 18446744078473872438), Vector[row](%id = 18446744078473872558), Vector[row](%id = 18446744078473872678), Vector[row](%id = 18446744078473872798), Vector[row](%id = 18446744078473872918), Vector[row](%id = 18446744078473873038), Vector[row](%id = 18446744078473873158), Vector[row](%id = 18446744078473873278), Vector[row](%id = 18446744078473873398), Vector[row](%id = 18446744078473873518), Vector[row](%id = 18446744078473873638), Vector[row](%id = 18446744078473874478), Vector[row](%id = 18446744078473874598), Vector[row](%id = 18446744078473874718), Vector[row](%id = 18446744078473874838), Vector[row](%id = 18446744078473874958), Vector[row](%id = 18446744078473875078), Vector[row](%id = 18446744078473875198), Vector[row](%id = 18446744078473875318), Vector[row](%id = 18446744078473867262), Vector[row](%id = 18446744078473867382)]

(3)

CurveBundle := proc()
  description "arg 1: fixed parameter\narg 2: fixed value\narg 3: color\narg 4: legend location";

  local minames, col, val, pts:

  minames := [`#mi("&epsilon;")`, `#mi("k")`];

  if _passed[1] = epsilon then
    col := 1:
  elif _passed[1] = k then
    col := 2:
  else
    error "first parameter passed should be epsilon or k"
  end if:

  if _passed[1] = epsilon
     and
     `not`(member(_passed[2], convert(epsilon_values, set))) then
     error cat(_passed[2], " is not a value in epsilon_values")
  end if:

  if _passed[1] = k
     and
     `not`(member(_passed[2], convert(k_values, set))) then
     error cat(_passed[2], " is not a value in k_values")
  end if:
  
  # same type of code could be written to ask for the value of epsilon or k
  # (I'm too lazzy to write it)

  val := _passed[2]:
  pts := map(ke -> if AllMax[ke, col]=val then convert(AllMax[ke, [3-col, 3]], list) end if, [$1..KE]);

  return plot(
           pts
           , color=_passed[3]
           , labels=[minames[3-col], "Max"]
           , legend=typeset(minames[col]=val)
           , legendstyle=[location=_passed[4]]
         )

end proc:
  

Describe(CurveBundle):


# arg 1: fixed parameter
# arg 2: fixed value
# arg 3: color
# arg 4: legend location
CurveBundle( )

 

CurveBundle(epsilon, 3.14, red, right);
PlotOneCurve(k, 3.14, red, right);

Error, (in CurveBundle) 3.14 is not a value in epsilon_values

 

Error, (in PlotOneCurve) 3.14 is not a value in k_values

 

RandColor := () -> ColorTools:-Color([seq(rand(0. .. 1.)(), i=1..3)]):

plots:-display(
  seq(CurveBundle(epsilon, val, RandColor(), right), val in epsilon_values)
  , size=[500, 400]
)

 

RandColor := () -> ColorTools:-Color([seq(rand(0. .. 1.)(), i=1..3)]):

plots:-display(
  seq(CurveBundle(k, val, RandColor(), right), val in k_values[[seq](1..numelems(k_values), 2)])
  , size=[600, 400]
)

 
 

 

Download Maximum_help_mmcdara.mw

@AHSAN 

Also need to plot the ratio of PP_Max to PP against x fro different values of ϵ by fiixing other parameters.

I see what PP is (an expression which depends only on one single indeterminate x)  but I can't see sime expression for PP_Max.

@mehdi jafari 

It is likely what @Alfred_F had in mind but I'm not sure

restart

eq := diff(y(x),x) = y(x) + exp(x)*exp(-3*x)/2 + int(exp(x+t)*y(t),t=0..x);
IC := y(0)=1;

diff(y(x), x) = y(x)+(1/2)*exp(x)*exp(-3*x)+int(exp(x+t)*y(t), t = 0 .. x)

 

y(0) = 1

(1)

deq := diff(eq, x)

diff(diff(y(x), x), x) = diff(y(x), x)-exp(x)*exp(-3*x)+int(exp(x+t)*y(t), t = 0 .. x)+exp(2*x)*y(x)

(2)

# Isolate the integral tem in deq:

rel := isolate(deq, select(has, indets(deq), int)[])

int(exp(x+t)*y(t), t = 0 .. x) = diff(diff(y(x), x), x)-(diff(y(x), x))+exp(x)*exp(-3*x)-exp(2*x)*y(x)

(3)

# Use 'rel' to rewrite eq without integral term

eq1 := eval(eq, rel);

diff(y(x), x) = y(x)+(3/2)*exp(x)*exp(-3*x)+diff(diff(y(x), x), x)-(diff(y(x), x))-exp(2*x)*y(x)

(4)

eq2 := simplify( isolate(eq1, diff(y(x), x$2)) )

diff(diff(y(x), x), x) = 2*(diff(y(x), x))-y(x)-(3/2)*exp(-2*x)+exp(2*x)*y(x)

(5)

# Provide an initial condition for D(y)(0) and 'dsolve'


 

Download idp.mw

 

upload your worksheet using the big green up arrow in the menu bar.

Without your code it is impossible to understand where the error comes from (please specify your Maple version too)

An example with Maple 2015  ExponentialFit.mw

Does your version enable computing directly the  R-squared and Adjusted R-squared statistics?

restart
with(Statistics):
X := Vector([1, 2, 3, 4, 5, 6], datatype=float):
Y := Vector([2, 3, 4.8, 10.2, 15.6, 30.9], datatype=float):
param := Fit(a+b*t+c*t^2, X, Y, t, output=parametervalues);

 

This works the same for  NonlinearFit.

@Kitonum 

Good, thanks

@janhardo 

I'm not a big fan neither (at least in the "duplicate question" case).

I'd found if better if the deleted question was deplaced into some temporary place and its original content replaced by an explicit warning to the OP telling it the reason why its question was removed, with a proposal to link it to another one.
The problem would be then to manage potential never-ending discussions between the OP and the moderator.

So the radicality of the present has some advantages.
All the more, that I don't remember of any OP who, after complaining about a deleted question and receiving an explanation as to why, moved heaven and earth to have it reintegrated in its original position.

@salim-barzani 

One of the most common reason why a question is deleted is that it is about a problem which is very closed to the one mentionned in an previous question asked a few hours before.

The decision of deleting the "second" question, or sometimes to transfering it into the thread relative to the "first" question, is taken by a moderator. And yes, I agree, it is a subjective decision (note that, as explained in the link below, it is possible to contest it). 
But from what I've seen since I've been on this site, this type of decision is often justified, even if being its victim is all but pleasant.
Note that you can post a new question titled "Why was my question deleted?" to get more explanations.

About Mapleprimes' policy see here Moderation

@Paras31 

The main point of my comment concerned the comparison of different numerical schemes and mentioned symplectic schemes as particularly suitable for the integration of Hamiltonian systems.
Your comment focuses only on the last (anecdotal) remark I made. However, I don't understand its content but, after all, that's not my business.

Good continuation 

@Al86 

Look at this ansatz_2_mmcdara.mw

In addition: here is a A_Toy_Problem.mw which explains how the perturbation method operates and how to perform it step-by-step in Maple.

Note: I used the package LargeExpressions in this filr, not because it is necessary (you coud use coeff/coeffs instead), but to help you understand how I proceeeded in file ansatz_2_mmcdara.mw.

It seems method=gear is the method that best conserves the total energy  (absolute error less than 10-9 over the range t=0..100, so about 103 times smaller than Rosenbrock)

Hamiltoninan_procedure_mmcdara.mw

About Gear's method:
In French Hairer   (Gear's method is named BDF in this paper)
In English Novik et al. (specificaly chapter V. PERFORMANCE OF THE ALGORITHMS)

More generaly:
If you are interested in energy conserving methods for hamitonian systems I advice you to give a look to this introductory Wikipedia  article.

At last: your ode contains the term k*diff(x(t), t) that you cancel out setting k=0.
I guess you are aware that the ode no longer defines a conservative system as soon as k <> 0...

I'm not sure the ansatz you use is is the one that should be used.

In the attached file, up to relation (ç), you will found how to get an expression of the form

sum(A[k]*epsilon^k, k=M..N)

where each A[k] represents the kth order differential equation.
With the ansatz you chosed M=-9 and N=11.

From relation (10) to relation (19) the solutions x[-1](t), x[0](t), x[1](t) are computed.
You will see they all have the sape formal expression p+q*exp(-t), and thus is  x(t) writes p+q*exp(-t) too.

It's hard to go further without having the initial conditions for x[-1](t), x[0](t), x[1](t), but  it seems the ansatz you chosed is not well suited for the original equation.
Indeed, assuming x(t) <> 0 for each t in the integration domain, the ode can be transform into what I named ode_2.
Collecting ode_2  wrt epsilon writes 

sum(B[k]*epsilon^k, k=0..4)

So why did you used an ansatz containing x[-1](t)/epsilon?

ansatz_1_mmcdara.mw

@C_R 

and thank you for the comment.
Doing this kind animation was quite funny and I did a few others in a row. The first one was quite time-consuming, but once you have understood how plottools (mainly) works it is not that difficult to do amazing things.

I completely agree with your last sentence, even though I'm no longer involved in using MapleSim.


A story:
For years, I've been trying to promote MapleSim in my company. Until it made a deal with a big actor (Dassult Systèmes) which provides turnkey solutions, or so they say. In particular  a product named Dymola (based on Modelica) that some colleagues wrongly considered as an alternative to MapleSim (they just saw the pictures and movies were quite close but forgot the CAS aspect behind).
Whatever, when it goes to big money MapleSoft doesn't hold a candle to aggressive distributors, specially when outsourcing becomes the rule.
In the end MapleSim was rejected and we are now dependent on Dassault's goodwill to develop, or not, what we really need. End of the story.

@Christian Wolinski 

First reference
Gauth  ( same problem and same result of yours [which is obviously correct])

Second reference
Doubtnut which gives answer 16
But this is not the same problem and @yangtheary made a mistake while reformulating the problem: the term power set makes all the difference!
(note the typo {x|x| < 3 ...} instead of  {x : |x| < 3 ...} on the screenshot below [the explanations the previous link provides prove t is indeed a typo]  }
 



 
In case @yangtheary would know how Maple can provide a solution, here is a way
 

restart

A := { isolve({x > -3, x < 3}) }

{{x = -2}, {x = -1}, {x = 0}, {x = 1}, {x = 2}}

(1)

B := A minus {{x = -1}}

{{x = -2}, {x = 0}, {x = 1}, {x = 2}}

(2)

 

R := map(t -> {(x, y)=(rhs(t[]), abs(rhs(t[])))}, B)

{{(x, y) = (-2, 2)}, {(x, y) = (0, 0)}, {(x, y) = (1, 1)}, {(x, y) = (2, 2)}}

(4)

P := combinat:-powerset(R):

numelems(P)

16

(5)

 


Download With_Maple.mw

 

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