mmcdara

@salim-barzani I'm sorry but th...

Answered: mmcdara 7891

December 16 2024

0 0

I'm sorry but the font is too small and the text too fuzzy for me to see clearly what it's written.
Is the Sawaqa-Kotera equation this one:

SK := diff(u, t) + 45*u^2*diff(u, x) - 15*diff(u, x)*diff(u, x$2) - 15*u*diff(u, x$3) + diff(u, x$5) = 0

Could you sent me something more readable?

Meanwhile here is the KdV equation with arbitrary alpha:

>

restart

PART ONE

>	alias(F=F(x, t), G=G(x, t))

(1)

>	with(PDEtools): undeclare(prime):

(2)

>	ND := proc(F, G, U) local v, w, f, g, a: v := op(F): if v[1] in U then w := -v[1] else w := v[1] end if: if v[2] in U then w := w, -v[2] else w := w, v[2] end if: f := op(0, F): g := op(0, G): a := diff(f(w)*g(v), U); convert(subs([w]=~[v], a), diff) end proc:

KdV with w = α ln (F)

PART TWO

A name like KdV__n refers to the equation (n) in the paper you present

>	alias(u=u(x, t), w=w(x, t))

(3)

>	KdV__1 := diff(u, x$3)+6udiff(u, x)+diff(u, t) = 0

(4)

>	KdV__3 := eval(KdV__1, u=diff(w, x$2))

(5)

>	KdV__4 := int~(lhs(KdV__3), x) = 0

(6)

>	prepa := eval(KdV__4, w=alpha*ln(F)): KdV__7 := op(-1, numer(lhs(prepa))) = 0;

3*alpha*(diff(diff(F, x), x))^2*F^2-6*alpha*(diff(diff(F, x), x))*(diff(F, x))^2*F+3*alpha*(diff(F, x))^4+(diff(diff(diff(diff(F, x), x), x), x))*F^3+(diff(diff(F, t), x))*F^3-4*(diff(diff(diff(F, x), x), x))*(diff(F, x))*F^2-(diff(F, t))*(diff(F, x))*F^2-3*(diff(diff(F, x), x))^2*F^2+12*(diff(diff(F, x), x))*(diff(F, x))^2*F-6*(diff(F, x))^4 = 0

(7)

PART THREE: express KdV__7 using ND "operators"

>	# As the highest derivation degree wrt x is 4, the D-operator will contain # ND(F, F, [x$4]): `#msubsup(mo("D"),mo("x"),mo("4"))` = ND(F, F, [x$4]);

(8)

>	ToRewrite := diff(F, x$4); RewriteAs := isolate((8), ToRewrite); Rewritten := simplify(subs(RewriteAs, lhs(KdV__7)));

diff(diff(diff(diff(F, x), x), x), x) = -(1/2)*(-8*(diff(diff(diff(F, x), x), x))*(diff(F, x))+6*(diff(diff(F, x), x))^2-`#msubsup(mo("D"),mo("x"),mo("4"))`)/F

3*alpha*(diff(diff(F, x), x))^2*F^2-6*alpha*(diff(diff(F, x), x))*(diff(F, x))^2*F+3*alpha*(diff(F, x))^4+(diff(diff(F, t), x))*F^3-(diff(F, t))*(diff(F, x))*F^2-6*(diff(diff(F, x), x))^2*F^2+12*(diff(diff(F, x), x))*(diff(F, x))^2*F-6*(diff(F, x))^4+(1/2)*F^2*`#msubsup(mo("D"),mo("x"),mo("4"))`

(9)

>	# Rewrite the terms containing the second derivative of F wrt x `#msubsup(mo("D"),mo("x"),mo("2"))` = ND(F, F, [x$2]);

(10)

>	ToRewrite := diff(F, [x$2]); RewriteAs := isolate((10), ToRewrite); Rewritten := simplify(subs(RewriteAs, Rewritten));

diff(diff(F, x), x) = -(1/2)*(-2*(diff(F, x))^2-`#msubsup(mo("D"),mo("x"),mo("2"))`)/F

(11)

>	# Rewrite the terms containing the second derivative of F wrt x and t `#msubsup(mo("D"),mrow(mo("x"),mo("t")))` = ND(F, F, [x, t]);

(12)

>	ToRewrite := diff(F, [x, t]); RewriteAs := isolate((12), ToRewrite); Rewritten := collect(simplify(expand(subs(RewriteAs, Rewritten))), [F, `#msubsup(mo("D"),mo("x"),mo("2"))`])

diff(diff(F, t), x) = -(1/2)*(-2*(diff(F, x))*(diff(F, t))-`#msubsup(mo("D"),mrow(mo("x"),mo("t")))`)/F

(13)

>	# The alpha=2 case simplify(eval(Rewritten, alpha=2)); # Operator form of KdV__7 Result(alpha=2) = ``(op(2,%))(F.F)

(14)

>	# Check: F^2 * ND(F, F, [x$4]) / 2 + F^2 * ND(F, F, [x, t]) / 2 + (3/4 * alpha - 3/2) * ND(F, F, [x$2])^2; simplify(% - lhs(KdV__7))

(15)

>	# Operator form of KdV__7 KdV__9 := ``(Rewritten)(F.F)=0

(16)

"To get a relation closer to equation (9) in your paper we need to prove this; D[xt] = D[x] @ D[t] = D[t] @ D[x]"

>

`#msub(mo("D"),mrow(mo("x"),mo("t")))` = ND(F, G, [x, t]);

# some algebraic transformations
     relation_x := `#msub(mo("D"),mo("x"))` = ND(F, G, [x]):
     alias(f=f(x, t), g=g(x, t)):
     equivalences := {diff(F, x)=f, diff(G, x)=g}:

     eval(relation_x, equivalences):
     map(z -> sign(z)*ND(op(sign(z)*z), [t]), [op(rhs(%))]):
     add(%):

`#msub(mo("D"),mo("x"))`@`#msub(mo("D"),mo("t"))` = eval(%, (rhs=lhs)~(equivalences));

# some algebraic transformations
     relation_t := `#msub(mo("D"),mo("t"))` = ND(F, G, [t]):
     equivalences := {diff(F, t)=f, diff(G, t)=g}:

     eval(relation_t, equivalences):
     map(z -> sign(z)*ND(op(sign(z)*z), [x]), [op(rhs(%))]):
     add(%):

`#msub(mo("D"),mo("t"))`@`#msub(mo("D"),mo("x"))` = eval(%, (rhs=lhs)~(equivalences))

(17)

Download Hirota_derivative_KdV_equation_2.mw

Formal derivation of equation (9)...

Answered: mmcdara 7891

December 16 2024

0 0

>

restart

PART ONE

>	alias(F=F(x, t), G=G(x, t))

(1)

>	with(PDEtools): undeclare(prime):

(2)

>	ND := proc(F, G, U) local v, w, f, g, a: v := op(F): if v[1] in U then w := -v[1] else w := v[1] end if: if v[2] in U then w := w, -v[2] else w := w, v[2] end if: f := op(0, F): g := op(0, G): a := diff(f(w)*g(v), U); convert(subs([w]=~[v], a), diff) end proc:

>	ND(F, G, [x]); ND(F, G, [t]);

(3)

>	ND(F, F, [x]); ND(F, F, [x, x]);

(4)

>	ND(F, G, [x$3]);

(5)

>	ND(F, G, [x, t]); `#msub(mo("D"),mrow(mo("x"),mo("t")))` = eval(%, G=F);

(6)

PART TWO

A name like KdV__n refers to the equation (n) in the paper you present

>	alias(u=u(x, t), w=w(x, t))

(7)

>	KdV__1 := diff(u, x$3)+6udiff(u, x)+diff(u, t) = 0

(8)

>	KdV__3 := eval(KdV__1, u=diff(w, x$2))

(9)

>	KdV__4 := int~(lhs(KdV__3), x) = 0

(10)

>	prepa := eval(KdV__4, w=2*ln(F)): KdV__7 := op(-1, numer(lhs(prepa))) = 0;

(11)

PART THREE: express KdV__7 using ND "operators"

>	# As the highest derivation degree wrt x is 4, the D-operator will contain # ND(F, F, [x$4]): `#msubsup(mo("D"),mo("x"),mo("4"))` = ND(F, F, [x$4]);

(12)

>	ToRewrite := diff(F, x$4)*F; RewriteAs := isolate((12), ToRewrite); Rewritten := subs(RewriteAs, lhs(KdV__7));

(13)

>	# As the highest total derivation degree is 2 and then ND(F, F, [x, t]) has to be involved `#msub(mo("D"),mrow(mo("x"),mo("t")))` = ND(F, F, [x, t]);

(14)

>	ToRewrite := F*(diff(F, [x,t])); RewriteAs := isolate((14), ToRewrite); Rewritten := subs(RewriteAs, Rewritten);

(15)

>	# Check: ND(F, F, [x$4])/2 + ND(F, F, [x, t])/2; simplify(% - lhs(KdV__7))

(16)

>	# Operator form of KdV__7 KdV__9 := ``(2*Rewritten)(F.F)=0

(17)

"To get a relation closer to equation (9) in your paper we need to prove this; D[xt] = D[x] @ D[t] = D[t] @ D[x]"

>

`#msub(mo("D"),mrow(mo("x"),mo("t")))` = ND(F, G, [x, t]);

# some algebraic transformations
     relation_x := `#msub(mo("D"),mo("x"))` = ND(F, G, [x]):
     alias(f=f(x, t), g=g(x, t)):
     equivalences := {diff(F, x)=f, diff(G, x)=g}:

     eval(relation_x, equivalences):
     map(z -> sign(z)*ND(op(sign(z)*z), [t]), [op(rhs(%))]):
     add(%):

`#msub(mo("D"),mo("x"))`@`#msub(mo("D"),mo("t"))` = eval(%, (rhs=lhs)~(equivalences));

# some algebraic transformations
     relation_t := `#msub(mo("D"),mo("t"))` = ND(F, G, [t]):
     equivalences := {diff(F, t)=f, diff(G, t)=g}:

     eval(relation_t, equivalences):
     map(z -> sign(z)*ND(op(sign(z)*z), [x]), [op(rhs(%))]):
     add(%):

`#msub(mo("D"),mo("t"))`@`#msub(mo("D"),mo("x"))` = eval(%, (rhs=lhs)~(equivalences))

(18)

>

Download Hirota_derivative_KdV_equation.mw

Gve more details......

Answered: mmcdara 7891

December 14 2024

0 0

@emrantohidikub

... but I can't guarantee that I'll be capable to answer your questions, or even that I'll have the time to do it if it requires a lot of work.

Rotational kinetic energy, contact, and ...

Answered: mmcdara 7891

December 14 2024

0 0

Not accounting for the rotational kinetic energy and assuming the ball has allway a single contact point with the surface is a huge approximation. Even in the absence of gravity (let's say the ball has an internal engine), a single contact requires that the curvature radius of the surface is everywhere larger than the radius of the ball.
More og that the displacement of the ball may exhibt pure rolling, pure sliding, or a a mix of these twos depending of the value of the contat force.

I aggree with Rouben saying it is quite a challenge to account for all these points in the most general way (even if there are probably a lot of papers online about this subject).

For information @C_R already did something close with MapleSim, see here

old fashion...

Answered: mmcdara 7891

December 13 2024

0 0

@acer

I asked years ago a question about the old fashion way Maple used to represent lengthy expression (using things like %1, %2, ... in the global expression and adding lines to the output of the form %1=something, %2=something-else, ...).
I couldn't find the answer I received but this could lead to something more readable here.

WHY...

Answered: mmcdara 7891

December 13 2024

0 0

didn't you just do what I suggested you HERE

MatlabFile := cat(currentdir(), "/STX.txt"):   # change the path and name as you want
ExportMatrix(MatlabFile, STX, target=MATLAB, format=rectangular, mode=ascii, format=entries);

The end of your code is

cat(currentdir()); # Which is a directory NOT a file
                    "C:\Program Files\Maple 2024"

# Below matlabData should be the name of a file: as its not you get an error
# The warning message results from STX being a matrix of complex.

ExportMatrix(matlabData, STX, target = MATLAB, format = rectangular, mode = ascii);

Warning, Matlab may not know how to read complex values stored in 'rectangular' format; try 'entries' format if problems arise
Error, (in ExportMatrix) permission denied

#  matlabData is unassigned
matlabData;
                 matlabData;

Begin to assign matlabData to a file with a legit absolute path.
Next use format=entries as I already told you (if you don't understand why read carefully the ExportMatrix help page): this will suppress the warning.

I advice you to read the help pages about cat and currentdir() too.

PS: Don't be surprised if a moderator deletes this question, which obviously duplicates the one I've already answered.

A suggestion...

Answered: mmcdara 7891

December 13 2024

0 0

@salim-barzani

Look at he ExportMatrix help page, items MATLAB® Binary Format or MATLAB® ASCII Format.
As you want to export a matrix of complex search (CTRL+F) for the word entries.

Here is an ASCII example

MatlabFile := cat(currentdir(), "/STX.txt"):   # change the path and name as you want
ExportMatrix(MatlabFile, STX, target=MATLAB, format=rectangular, mode=ascii, format=entries);

As I don't have Matlab at home I can't check it can read and plot the exported matrix correctly.

In case you want avoid using the Finance package, here is a Wiener process generator

>

restart

>	with(Statistics):

>	G := RandomVariable(Normal(0, 1))

(1)

>	N := 1000: T := 10: dt := T/N: S := Sample(G, N): dW := sqrt(dt) *~ S: W := Array(<0>): C := CumulativeSum(dW): W := ArrayTools:-Extend(W, CumulativeSum(dW), inplace=true): plot([seq(0..T, dt)], W, style=line, size=[1000, 400]);

>

# Several Wiener processes

rep := 20:
all := NULL:
col := () -> ColorTools:-Color([seq(rand(0. .. 1.)(), i=1..3)]):

for i from 1 to rep do
  S := Sample(G, N):
  dW := sqrt(dt) *~ S:
  W := Array(<0>):
  C := CumulativeSum(dW):
  W := ArrayTools:-Extend(W, CumulativeSum(dW), inplace=true):

  all := all, plot([seq(0..T, dt)], W, style=line, size=[1000, 400], color=col()):
end do:

plots:-display(all)

>

Download WienerProcessesGenerator.mw

Your initial code with W(t) plugged in...

Answered: mmcdara 7891

December 12 2024

0 0

@salim-barzani

Is this what you want (left = Re, right = Im)
graph-stochastic_mmcdara.mw

@dharr Thanks for your explanation....

Answered: mmcdara 7891

December 06 2024

0 0

@dharr

Thanks for your explanation.
See you next time.

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

@Ronan The denomintor then is&nbs...

Answered: mmcdara 7891

December 06 2024

0 0

@Ronan

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

For smart display only...

Answered: mmcdara 7891

November 25 2024

0 0

>

restart:

>	edp := diff(u(x, t), t)+lambda*diff(u(x, t), x$2) = 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);

(2)

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

(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);

(4)

Download smart_display.mw

Data file...

Answered: mmcdara 7891

November 25 2024

0 0

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

Answered: mmcdara 7891

November 23 2024

0 0

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

Answered: mmcdara 7891

November 19 2024

0 0

@segfault

Good, I'm glad I could help.

Can you pay attention to what you write?...

Answered: mmcdara 7891

November 19 2024

0 0

@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

(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];

(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]

(3)

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

(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")
);

>

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")
);

>

Download Accept_it_or_not_but_I_am_done.mw

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post

7891 Reputation

22 Badges

MaplePrimes Activity

These are replies submitted by mmcdara

@salim-barzani I'm sorry but th...

Formal derivation of equation (9)...

Gve more details......

Rotational kinetic energy, contact, and ...

old fashion...

WHY...

A suggestion...

Your initial code with W(t) plugged in...

@dharr Thanks for your explanation....

@Ronan The denomintor then is&nbs...

For smart display only...

Data file...

@Carl Love @dharr @acerT...

@segfault Good, I'm glad I coul...

Can you pay attention to what you write?...

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Save this setting as your default sorting preference?
Note: You can change your preference any time in your account settings
Don't show this again