mmcdara

6740 Reputation

18 Badges

8 years, 188 days

MaplePrimes Activity


These are replies submitted by mmcdara

@dharr 

Thanks for your explanation.
See you next time.

By the way: congratulations, you've passed me :-)

@Ronan 

The denomintor then is 100*n^(33/20).

restart:

edp := diff(u(x, t), t)+lambda*diff(u(x, t), x$2) = s

diff(u(x, t), t)+lambda*(diff(diff(u(x, t), x), x)) = s

(1)

eDp := convert(edp, D):
inD := select(has, indets(eDp), D):
rw  := map(i -> i =  eval(op([0, 0], i), {D=u, 1=x, 2=t}), inD):

rwedp := eval(eDp, rw);

lambda*u[x, x]+u[t] = s

(2)

original_edp := convert(eval(rwedp, (rhs=lhs)~(rw)), diff)

diff(u(x, t), t)+lambda*(diff(diff(u(x, t), x), x)) = s

(3)

# Or, if you want to get rid of the commas:

rw  := map(i -> i = cat('u__', eval(op(op([0, 0], i)), {1=x, 2=t})), inD):
rwedp := eval(eDp, rw);

lambda*u__xx+u__t = s

(4)
 

 

Download smart_display.mw

Is located in your own directories.
Please upload the file or provide us some link we can acccess.

Alread: writting this in your procedure SSE is illegal

.
.
.
local solution := dsolve({ode1,ode2,S(0) = S0,T(0) = i0},parameters=[beta,Gamma,S0,i0],[S(t),T(t)],numeric);
solution(parameters = [b,g,x,y]);
local difference := add((rhs(solution(k)[3])-sum_total[k])^2,k=1..N);
return difference;
end proc:

A correct syntaxt is

.
.
.
local solution := dsolve({ode1,ode2,S(0) = S0,T(0) = i0},parameters=[beta,Gamma,S0,i0],[S(t),T(t)],numeric);
solution(parameters = [b,g,x,y]);
add((rhs(solution(k)[3])-sum_total[k])^2,k=1..N);
end proc:


Last but not least: if you expect finding 4 parameters such that the brown curve is close to the empirical data, this is a dead end: don't need to be a great scientist to understand why.

@Carl Love  @dharr @acer

Thanks you all for your involvement.
I particularly appreciated Carl's explanations and embedded comments by Dharr (I missed the usedpts=1044 for method=_MonteCarlo).

There are stlll a few messages that bother me.

For instance method=_CubaSuave returns 
cuba: chi-square probability that the error is not reliable: 1... which means that there are 100% of chances that the error (cuba: estimated (absolute) error: 6.73167e-05) ! Are you sure the first message is correct?

All the more that method=_CubaVegas gives
cuba: estimated (absolute) error: 6.58857e-05
cuba: chi-square probability that the error is not reliable: 0.0023283

and that CubaSuave is aimed to be an improvement of CubaVegas (see  Hahn for instance)

By the way, as CubaVegas uses a Quasi Random Numbers Generator [QRNG] (the Sobol' one, see previous ref) don't you think it would have been a good idea for the core development team to also offer this quasi-random number generator (see NAG library QMC) ?

(For the record there exist several QRNG and I implemented Faure's in Maple a few years ago. But I believe fast built-in functions for QRNG would be welcome.)

 

@segfault 

Good, I'm glad I could help.

@Aung 

YourModel :

epsilon[0] + epsilon[0] * alpha[theta] * B[1] * (-1 + exp(-beta[1] * t)) 
           + epsilon[0] * alpha[theta] * B[2] * (-1 + exp(-beta[2] * t))
           + epsilon[0] * alpha[theta] * B[3] * (-1 + exp(-beta[3] * t))

Where  epsilon[0] > 0 and alpha[theta] > 0
As you impose B[1] >= 0, B[2] >= 0, B[3] >= 0, beta[1] >= 0, beta[2] >= 0, beta[3] >= 0, 
the 3 last terms are always negative (t >=0) or null in the "best" case.
Given that epsilon[0] = 0.5581381911e-2 the maximum value of your model cannot exceed this value:
So how can you even think this model could represent your data? 




My advice, and do with it what you will, is summarized in the attached file:
(the model I use is similar to yours after sign correction)
 

restart;

with(Statistics):with(plots):with(Optimization):with(LinearAlgebra):


# given data from strain rate curve
E_0[theta] := 7.883352314*10^9;
alpha[theta]:= 0.982

7883352314.

 

.982

(1)


# experimental creep data under 44 at 100 degree celcius
c_strain := Vector ([<<0>,<0.0284698>,<0.0533808>,<0.0782918>,<0.0996441>,<0.124555>,<0.142349>,<0.156584>,<0.16726>,<0.177936>,<0.181495>,<0.188612>,<0.192171>,<0.19573>,<0.19573>,<0.202847>,<0.206406>,<0.206406>,<0.209964>,<0.209964>,<0.209964>,<0.206406>,<0.209964>>]):

c_time := Vector ([<<0>,<0>,<0.048>,<0.192>,<0.352>,<0.544>,<0.704>,<0.896>,<1.088>,<1.312>,<1.52>,<1.76>,<1.984>,<2.208>,<2.464>,<2.736>,<3.088>,<3.392>,<3.664>,<4.016>,<4.352>,<4.592>,<4.832>>]):
sigma[0] := 44*10^6;
epsilon[0] := sigma[0]/E_0[theta];

44000000

 

0.5581381911e-2

(2)


# change vector to list
c_strain := convert(c_strain,list):
c_time := convert(c_time,list):

# extract zero from list
c_strain := c_strain [2..-1];
c_time := c_time [2..-1];

[0.284698e-1, 0.533808e-1, 0.782918e-1, 0.996441e-1, .124555, .142349, .156584, .16726, .177936, .181495, .188612, .192171, .19573, .19573, .202847, .206406, .206406, .209964, .209964, .209964, .206406, .209964]

 

[0, 0.48e-1, .192, .352, .544, .704, .896, 1.088, 1.312, 1.52, 1.76, 1.984, 2.208, 2.464, 2.736, 3.088, 3.392, 3.664, 4.016, 4.352, 4.592, 4.832]

(3)


# for further calculation need to know how many elements are in the list
M := nops(c_strain);
N := nops(c_time);

22

 

22

(4)


WHAT I SUGGESTED YOU

CreepCompliance := (t, b) ->  (1 -exp(-b*t)) / b;

NComponentModel := n -> unapply( shift + add((s||i)*CreepCompliance(t, b||i) , i=1..n), t):

NComponentModel(3);

eps := 1e-6:

Phi := 0;
obj := add( ( NComponentModel(3)~(c_time) - c_strain )^~2 )
           +
           Phi / (1e-4+add( (s||i - add(s||j, j=1..3)/3)^2, i=1..3)):

 opt_3 := Optimization:-NLPSolve( obj, {seq(s||i >= eps, i=1..3), seq(b||i >= eps, i=1..3)}, iterationlimit=1000);

model_3 := eval(NComponentModel(3)(t), (opt_3)[2]):

[ seq(isolate(epsilon[0]*alpha[theta]*B[i] = eval(s||i / b||i, opt_3[2]), B[i]), i=1..3) ];

plots:-display(
   Statistics:-ScatterPlot(c_time, c_strain, symbol=circle)
   , plot(model_3, t=0. .. max(c_time), color=red,  legend="3 components")
);

proc (t, b) options operator, arrow; (1-exp(-b*t))/b end proc

 

proc (t) options operator, arrow; shift+s1*(1-exp(-b1*t))/b1+s2*(1-exp(-b2*t))/b2+s3*(1-exp(-b3*t))/b3 end proc

 

0

 

[0.225314772163480432e-3, [b1 = HFloat(1.3084255135648692), b2 = HFloat(1.3084255135598823), b3 = HFloat(1.308425513559216), s1 = HFloat(0.07458889456019333), s2 = HFloat(0.07458889455231471), s3 = HFloat(0.07458889455084408), shift = HFloat(0.03700542692440384)]]

 

[B[1] = HFloat(10.400924369210232), B[2] = HFloat(10.400924368151253), B[3] = HFloat(10.40092436795148)]

 

 

Phi := 1;
obj := add( ( NComponentModel(3)~(c_time) - c_strain )^~2 )
           +
           Phi / (1e-4+add( (s||i - add(s||j, j=1..3)/3)^2, i=1..3)):

 opt_3 := Optimization:-NLPSolve( obj, {seq(s||i >= eps, i=1..3), seq(b||i >= eps, i=1..3)}, iterationlimit=1000);

model_3 := eval(NComponentModel(3)(t), (opt_3)[2]):

[ seq(isolate(epsilon[0]*alpha[theta]*B[i] = eval(s||i / b||i, opt_3[2]), B[i]), i=1..3) ];

plots:-display(
   Statistics:-ScatterPlot(c_time, c_strain, symbol=circle)
   , plot(model_3, t=0. .. max(c_time), color=red,  legend="3 components")
);

1

 

[0.934443784351978465e-4, [b1 = HFloat(370724.4386808285), b2 = HFloat(1.2325972978366215), b3 = HFloat(817606.2552664457), s1 = HFloat(4721.29039619488), s2 = HFloat(0.2034032484701095), s3 = HFloat(2179.004303789217), shift = HFloat(0.028462609697833815)]]

 

[B[1] = HFloat(2.3235727208121157), B[2] = HFloat(30.108106594130508), B[3] = HFloat(0.48625116430224813)]

 

 

 

 


 

Download Accept_it_or_not_but_I_am_done.mw

@C_R 

Let F =~ 0.
Convert F into exponential and apply the substitution

X2Y := { exp(I*x4) = y4, exp((I*x5) = y5}

to the result.
This will avoid the square roots.
Then F is a rational fraction whose denominator, equal to y4*y5, is not null (y4 and Y( are complex numbers). So consider the numerator N (N =~ 0)
Maybe this could this could make Isolate's task easier ???

Another idea: instead of solving {F1, F2, F3, F4, F5} try solving {F1, F2-F1, F3-F1, F4-F1, F5-F1}

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

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