mmcdara

output=Array...

Answered: mmcdara 7072

December 24 2024

1 0

The simplest way, IMO, to get this table is to use the option output=Array in dsolve/numeric:

when    := [seq(0..50, 1)]:
sol     := dsolve(`union`(sys, ic), numeric, output=Array(when)):
sol[1];  # to verify the order of the outputs
results := round~(sol[2][1]):

printf("%-10s %-15s %-15s %-15s\n", "Day", "Infected", "Recovered", "Susceptible");
printf("---------------------------------------------\n");

fmt := "%-10d %-15d %-15d %-15d\n":
map(i -> printf(fmt, entries(results[i+1], nolist)), when):

Another way...

Answered: mmcdara 7072

December 23 2024

4 1

y := substring(x, 1..1)

To get the value of y:

x := asdf:
a := 3:
y := substring(x, 1..1);
                               a
eval(y)
                               3

For information only...

Answered: mmcdara 7072

December 22 2024

0 2

This doesn't answer your question, just to say that geom3d provides an elegant way (IMO) to determine if a point is on a line.

restart
with(geom3d):
line(L, [3+2*alpha, 1+6*alpha, 4-5*alpha], alpha):

point(A, 9, 19, -11):
IsOnObject(A, L)
                              true
point(B, 9, 1, -11):
IsOnObject(B, L)
                             false

Use similar triangles...

Answered: mmcdara 7072

December 22 2024

3 0

Plus an animation to illustrate the demonstration
Similar_Triangles_Animation.mw

A less brute force...

Answered: mmcdara 7072

December 20 2024

2 1

Here is a much faster (still quite brute) way to get the result

>

restart:

>	N := n -> add(n[i]*10^(i-1), i=1..6):

>	P := [ seq(ListTools:-Rotate([seq(n[i], i=1..6)], j), j=0..5) ]

[[n[1], n[2], n[3], n[4], n[5], n[6]], [n[2], n[3], n[4], n[5], n[6], n[1]], [n[3], n[4], n[5], n[6], n[1], n[2]], [n[4], n[5], n[6], n[1], n[2], n[3]], [n[5], n[6], n[1], n[2], n[3], n[4]], [n[6], n[1], n[2], n[3], n[4], n[5]]]

(1)

>	tstart := time(): NP := N~(P): Z := NP[2..-1]: ZP := combinat:-permute(Z): sols := NULL: for m from 1 to numelems(ZP) do eval(isolve({seq(NP[1]*(k+1) = ZP[m][k], k=1..5)}), _Z1=1); if max(rhs~(%)) < 10 then sols := sols, % end if; end do: sols; time()-tstart;

(2)

>	N([eval(seq(n[i], i=1..6), sols)]); % *~ [$1..6];

(3)

>	# vv code tstart := time(): P:=combinat:-permute([seq](0..9),6): # nops(%);

>	for L in P do

>	n:=parse(cat(L[]));

>

ok:=true;

>	for k from 2 to 6 do

>

n1:=n*k;

>	L1:=convert(n1,base,10):

>	if {L[]} <> {L1[]} then ok:=false; break fi;

>

od;

>	if ok then print(n, x23456=[seq(k*n,k=2..6)]) fi;

>	od: time()-tstart

>

(4)

Download Improvement.mw

Here it is...

Answered: mmcdara 7072

December 20 2024

2 5

After a lot of work as you didn't provide all the necessary details:

Note there is likely an error in the paper you present: the transformation should probably be

u = -2 * ln(f_xx) instead of u = 2 * ln(f_xx)

, as this screen capture seems to confirm

>

restart

Emulation of the Hirota derivation operator

>	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:

Verify if
ND(F, F, [x, t]) + ND(F, F, [x$6])
is indeed equal to the lhs of equation (3.2) in the excerpt of the paper you present.

>	ND(f, f, [x, t]) + ND(f, f, [x$6]);

(3)

Sawada-Kotera equation

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

(4)

>	SK__u := diff(u, t) + 45u^2diff(u, x) - 15diff(u, x)diff(u, x$2) - 15udiff(u, x$3) + diff(u, x$5) = 0

(5)

>	SK__f := eval(SK__u, u=diff(f, x$2))

(6)

>	ISK__f := map(Int, lhs(SK__f), x) = 0

(7)

>	ISK__fv := value(ISK__f);

(8)

The problem here is the uneval integral term.
Let Z its integrand:

>	Z := IntegrationTools:-GetIntegrand( select(has, [op(lhs(ISK__fv))], int)[] );

(9)

Use two by part integrations to compute its integral:

>	use IntegrationTools in ``(Int(Z, x)) = Parts(Int(Z, x), diff(f, x$2)); lhs(%) = expand(op(1, rhs(%)) + Parts(Expand(op(2, rhs(%))), diff(f, x$3))); % + %%; map(Expand, %); Reductor := expand~(% /~ 2); end use;

(10)

Use this relation to reduce ISK__f and get an expression free from integrals:

>	eval(IntegrationTools:-Expand(ISK__f), Reductor); ISK__fv := value(%);

(11)

Use the transformation F = -2*ln(f)

>	alias(F = F(x, t)): ISK__F := normal(eval(ISK__fv, f=alpha*ln(F))): ISK__F := eval(numer(lhs(ISK__F)), alpha=-2) = 0;

(12)

>	# As the highest derivation degree wrt x is 6, the D-operator will contain # ND(F, F, [x$6]): x6 := `#msubsup(mo("D"),mo("x"),mo("6"))` = ND(F, F, [x$6]); ToRewrite := diff(F, x$6): RewriteAs := isolate(x6, ToRewrite); Rewritten := simplify(subs( eval(RewriteAs, alpha=2), lhs(ISK__F)));

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

(13)

>

# Rewrite the terms containing the second derivative of F wrt x and t

xt := `#msubsup(mo("D"),mrow(mo("x"),mo("t")))` = ND(F, F, [x, t]);

ToRewrite := diff(F, [x, t]):
RewriteAs := isolate(xt, 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, t))*(diff(F, x))-`#msubsup(mo("D"),mrow(mo("x"),mo("t")))`)/F

(14)

>	SK_H := ``(eval(-Rewritten, F=1))(F.F)=0

(15)

"To get a relation closer to the las one in the excerpt of the paper you present 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))

(16)

Then

>	SK_H := eval((15), { lhs(xt) = `#msub(mo("D"),mo("x"))` * `#msub(mo("D"),mo("t"))`, F=f })

(17)

Download Sawada_Kotera_(Hirota_form).mw

'known' option...

Answered: mmcdara 7072

December 19 2024

3 3

@rcorless already answered you.

Nevertheless, if for some personal reason you prefer to solve a first system with unknown functions x(t), y(t) and z(t), and then plug the solutions within a second system with unknown functions P(t), Q(t) and R(t), you can use the known option of dsolve/numeric (look to the help page).
This is a very powerful feeature which is worth being known as it can be use in several other situations.

Here is how to use it

>

restart:

>	eq1 := diff(x(t), t)-(1/6)(6x(t)^3y(t)+(2y(t)^2-2)x(t)^2+3y(t)(z(t)-2)x(t)-2y(t)^2+2)sqrt(3) = 0;

>

>	eq2 := diff(y(t), t)-(1/6)(y(t)-1)sqrt(3)(y(t)+1)(6x(t)^2+2y(t)x(t)+3z(t)-2) = 0;

>

>	eq3 := diff(z(t), t)-(1/3)z(t)sqrt(3)(6y(t)x(t)^2+2x(t)y(t)^2+3z(t)y(t)-2x(t)-3*y(t)) = 0;

>

>	sys := eq1, eq2, eq3:

>	ics := x(0) = -0.01, y(0) = .99, z(0) = 0.01;

>

>	sol := dsolve({ics, sys}, type = numeric);

(1)

>	eq4 := diff(R(t), t)-P(t)Z(t)-(-2(-Y(t)^2+2)X(t)/sqrt(3)+sqrt(3)(-2X(t)^2+Z(t)+4/3)Y(t))*R(t) = 0;

>	eq5 := diff(Q(t), t)-(2/3)R(t)+2((1/3)Y(t)+X(t))P(t)/sqrt(3)-(-2(-Y(t)^2+2)X(t)/sqrt(3)+2sqrt(3)(X(t)^2-(1/2)Z(t)-2/3)X(t))*Q(t) = 0;

>	eq6 := diff(P(t), t)+(1/2)R(t)+2sqrt(3)X(t)Q(t)+(2(-Y(t)^2+2)X(t)/sqrt(3)+sqrt(3)(-2X(t)^2+Z(t)+1)Y(t))P(t) = 0;

>

>	SYS := eq4, eq5, eq6:

diff(Q(t), t)-(2/3)*R(t)+(2/3)*((1/3)*Y(t)+X(t))*P(t)*3^(1/2)-(-(2/3)*(-Y(t)^2+2)*X(t)*3^(1/2)+2*3^(1/2)*(X(t)^2-(1/2)*Z(t)-2/3)*X(t))*Q(t) = 0

diff(P(t), t)+(1/2)*R(t)+2*3^(1/2)*X(t)*Q(t)+((2/3)*(-Y(t)^2+2)*X(t)*3^(1/2)+3^(1/2)*(-2*X(t)^2+Z(t)+1)*Y(t))*P(t) = 0

(2)

>	indets({SYS});

(3)

>	# You didn't specicify ics ==> I use arbitrary ones ICS := P(0) = rand(0. .. 1.)(), Q(0) = rand(0. .. 1.)(), R(0) = rand(0. .. 1.)()

(4)

>

X := proc(s)
  if s::numeric then
    eval(x(t), sol(s))
  else
    'procname'(s)
end if:
end proc:

Y := proc(s)
  if s::numeric then
    eval(y(t), sol(s))
  else
    'procname'(s)
end if:
end proc:

Z := proc(s)
  if s::numeric then
    eval(z(t), sol(s))
  else
    'procname'(s)
end if:
end proc:

>	SOL := dsolve({SYS, ICS}, type=numeric, known=[X(t), Y(t), Z(t)])

(5)

>	# All in a row soln := dsolve({eq1, eq2, eq3, eq4, eq5, eq6, ics, ICS}, {P(t), Q(t), R(t), x(t), y(t), z(t)}, numeric, start = 0, range = 0 .. 1)

(6)

>	plots:-display( plots:-odeplot(SOL, [t, P(t)], t=0..1, color=red, thickness=3, legend="two systems"), plots:-odeplot(soln, [t, P(t)], t=0..1, color=blue, thickness=7, transparency=0.8, legend="one system") )

Download known_option.mw

Do you mean something like this?...

Answered: mmcdara 7072

December 16 2024

1 2

(sorry the content cannot be loaded, likely due to some plot command)

Maybe_this.mw

Illustration : epsilon_r for N=1..10 and P from 0.1 to 10

Nice (?) splitting of a latex expression...

Answered: mmcdara 7072

December 15 2024

0 4

@salim-barzani

As an example I used your lengthy eq3 expression without any simplifications, just to work on a big expression.

The idea is to split the Maple expression on several lines by adding break lines between operators in order ot introduce a smart breaking.
Your lengthy expression is then made of several sub-expressions.
It suffices now to generate the LaTeX code for each sub-expression, to concatenate to of those code the "\\\\" breaking line command, finally to concatenate all those codes.

NOTE: the `#mo("1",mathcolor="white")` trick enables having a subexpression beginning with a '+' character (which would be automaticlly removed by Maple)

>

restart

>

eq3 := (4*(-(2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda))*lambda+4*mu^2))*alpha[1]^2*a[5]*alpha[0]-12*mu^2*alpha[1]^2*a[5]*alpha[0]+(3*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda))*alpha[1]^2*alpha[0]*a[2]-(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*k^2*a[1]*alpha[1]^2+(1/2)*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*lambda*a[1]+(5*(-(2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda))*lambda+4*mu^2))*alpha[1]^4*alpha[0]*a[4]+(10*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda))*alpha[1]^2*alpha[0]^3*a[4]+(6*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda))*alpha[1]^2*alpha[0]^2*a[3]-6*lambda*beta[0]^2*alpha[1]^2*a[3]-2*lambda*beta[0]^2*a[5]*alpha[0]+6*mu*beta[0]*alpha[1]^2*a[2]+3*mu*beta[0]*a[5]*alpha[0]^2-9*mu^2*alpha[1]^2*a[1]*(1/4)+3*mu*a[1]*alpha[0]*beta[0]*(1/2)+(4*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda))*alpha[1]^2*lambda*a[5]*alpha[0]-20*mu*beta[0]*lambda*alpha[1]^4*a[4]+24*mu*beta[0]*alpha[1]^2*alpha[0]*a[3]-30*lambda*beta[0]^2*alpha[1]^2*alpha[0]*a[4]+60*mu*beta[0]*alpha[1]^2*alpha[0]^2*a[4]-7*mu*beta[0]*lambda*a[5]*alpha[1]^2-w*beta[0]^2-(1/4)*lambda*beta[0]^2*a[1]+3*beta[0]^2*alpha[0]*a[2]-k^2*a[1]*beta[0]^2+10*beta[0]^2*alpha[0]^3*a[4]+6*beta[0]^2*alpha[0]^2*a[3]-(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*w*alpha[1]^2+(1/4)*(3*(-(2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda))*lambda+4*mu^2))*alpha[1]^2*a[1]+(-(2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda))*lambda+4*mu^2)*alpha[1]^4*a[3] = 0;

(1)

>

LENGTH := 330: # to adjust to the page width in your latex file

printlevel:=2:
P := NULL:
L := [op(lhs(eq3))]:
k := 0:
while L <> [] do
  k := k+1:
  l := 0:
  i := 1:
  while length(l) < LENGTH and i <= numelems(L) do
    l := l + L[i]:
    i := i+1:
  end do:
  P := P, l:
  L := L[i..-1]:
end do:

P := [P]:
for i from 1 to numelems(P) do
  p := P[i]:
  if substring(convert(p, string), 1..1) <> "-" then
    P[i] := `#mo("1",mathcolor="white")` + p
  end if:
end do:

P[-1] := P[-1] = 0:
print~(P):

-(9/4)*mu^2*alpha[1]^2*a[1]+3*beta[0]^2*alpha[0]*a[2]+10*beta[0]^2*alpha[0]^3*a[4]+6*beta[0]^2*alpha[0]^2*a[3]-(1/4)*lambda*beta[0]^2*a[1]+(1/4)*(-6*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+12*mu^2)*alpha[1]^2*a[1]-w*beta[0]^2

`#mo("1",mathcolor="white")`+4*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*lambda*a[5]*alpha[0]-20*mu*beta[0]*lambda*alpha[1]^4*a[4]+24*mu*beta[0]*alpha[1]^2*alpha[0]*a[3]-30*lambda*beta[0]^2*alpha[1]^2*alpha[0]*a[4]+60*mu*beta[0]*alpha[1]^2*alpha[0]^2*a[4]

-7*mu*beta[0]*lambda*a[5]*alpha[1]^2-(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*w*alpha[1]^2+(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^4*a[3]+4*(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^2*a[5]*alpha[0]

-12*mu^2*alpha[1]^2*a[5]*alpha[0]+3*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]*a[2]+(1/2)*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*lambda*a[1]+5*(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^4*alpha[0]*a[4]

`#mo("1",mathcolor="white")`+10*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]^3*a[4]+6*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]^2*a[3]-6*lambda*beta[0]^2*alpha[1]^2*a[3]-2*lambda*beta[0]^2*a[5]*alpha[0]+6*mu*beta[0]*alpha[1]^2*a[2]

`#mo("1",mathcolor="white")`+3*mu*beta[0]*a[5]*alpha[0]^2+(3/2)*mu*a[1]*alpha[0]*beta[0]-36*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*a[1]*alpha[1]^2-36*a[1]*beta[0]^2 = 0

(2)

>

LENGTH := 150: # to adjust to the page width in your latex file

printlevel:=2:
P := NULL:
L := [op(lhs(eq3))]:
k := 0:
while L <> [] and k < 20 do
  k := k+1:
  l := 0:
  i := 1:
  while length(l) < LENGTH and i <= numelems(L) do
    l := l + L[i]:
    i := i+1:
  end do:
  P := P, l:
  L := L[i..-1]
end do:

P := [P]:
for i from 1 to numelems(P) do
  p := P[i]:
  if substring(convert(p, string), 1..1) <> "-" then
    P[i] := `#mo("1",mathcolor="white")` + p
  end if:
end do:

P[-1] := P[-1] = 0:
print~(P):

-(1/4)*lambda*beta[0]^2*a[1]+(1/4)*(-6*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+12*mu^2)*alpha[1]^2*a[1]

-w*beta[0]^2+4*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*lambda*a[5]*alpha[0]-20*mu*beta[0]*lambda*alpha[1]^4*a[4]

-7*mu*beta[0]*lambda*a[5]*alpha[1]^2-(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*w*alpha[1]^2+(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^4*a[3]

`#mo("1",mathcolor="white")`+4*(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^2*a[5]*alpha[0]-12*mu^2*alpha[1]^2*a[5]*alpha[0]

`#mo("1",mathcolor="white")`+3*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]*a[2]+(1/2)*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*lambda*a[1]

`#mo("1",mathcolor="white")`+5*(-2*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*lambda+4*mu^2)*alpha[1]^4*alpha[0]*a[4]+10*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]^3*a[4]

`#mo("1",mathcolor="white")`+6*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*alpha[1]^2*alpha[0]^2*a[3]-6*lambda*beta[0]^2*alpha[1]^2*a[3]

-36*(lambda*B[1]^2-lambda*B[2]^2-mu^2/lambda)*a[1]*alpha[1]^2-36*a[1]*beta[0]^2 = 0

(3)

>	# LaTeX output cat(map(p -> cat(latex(p, output=string), "\\\\"), P)[]); # If you have to remove the \mbox {{\tt `\#mo(\"1\",mathcolor=\"white\")`}} text: StringTools:-SubstituteAll(%, "\\mbox {{\\tt `\\#mo(""1"",mathcolor=""white"")`}}", ""):

$"-9/4\,{\mu}^{2}{\alpha_{{1}}}^{2}a_{{1}}+3\,{\beta_{{0}}}^{2}\alpha_{{0}}a_{{2}}+10\,{\beta_{{0}}}^{2}{\alpha_{{0}}}^{3}a_{{4}}+6\,{\beta_{{0}}}^{2}{\alpha_{{0}}}^{2}a_{{3}}\\-1/4\,\lambda\,{\beta_{{0}}}^{2}a_{{1}}+1/4\, \left( -6\, \left( \lambda\,{B_{{1}}}^{2}-\lambda\,{B_{{2}}}^{2}-{\frac {{\mu}^{2}}{\lambda}} \right) \lambda+12\,{\mu}^{2} \right) {\alpha_{{1}}}^{2}a_{{1}}\\-w{\beta_{{0}}}^{2}+4\, \left( \lambda\,{B_{{1}}}^{2}-\lambda\,{B_{{2}}}^{2}-{\frac {{\mu}^{2}}{\lambda}} \right) {\alpha_{{1}}}^{2}\lambda\,a_{{5}}\alpha_{{0}}-20\,\mu\,\beta_{{0}}\lambda\,{\alpha_{{1}}}^{4}a_{{4}}\\-30\,\lambda\,{\beta_{{0}}}^{2}{\alpha_{{1}}}^{2}\alpha_{{0}}a_{{4}}+60\,\mu\,\beta_{{0}}{\alpha_{{1}}}^{2}{\alpha_{{0}}}^{2}a_{{4}}+24\,\mu\,\beta_{{0}}{\alpha_{{1}}}^{2}\alpha_{{0}}a_{{3}}\\-7\,\mu\,\beta_{{0}}\lambda\,a_{{5}}{\alpha_{{1}}}^{2}- \left( \lambda\,{B_{{1}}}^{2}-\lambda\,{B_{{2}}}^{2}-{\frac {{\mu}^{2}}{\lambda}} \right) w{\alpha_{{1}}}^{2}+ \left( -2\, \left( \lambda\,{B_{{1}}}^{2}-\lambda\,{B_{{2}}}^{2}-{\frac {{\mu}^{2}}{\lambda}} \right) \lambda+4\,{\mu}^{2} \right) {\alpha_{{1}}}^{4}a_{{3}}\\\mbox {{\tt `\#mo("1",mathcolor="white")`}}+4\, \left( -2\, \left( \lambda\,{B_{{1}}}^{2}-\lambda\,{B_{{2}}}^{2}-{\frac {{\mu}^{2}}{\lambda}} \right) \lambda+4\,{\mu}^{2} \right) {\alpha_{{1}}}^{2}a_{{5}}\alpha_{{0}}-12\,{\mu}^{2}{\alpha_{{1}}}^{2}a_{{5}}\alpha_{{0}}\\\mbox {{\tt `\#mo("1",mathcolor="white")`}}+3\, \left( \lambda\,{B_{{1}}}^{2}-\lambda\,{B_{{2}}}^{2}-{\frac {{\mu}^{2}}{\lambda}} \right) {\alpha_{{1}}}^{2}\alpha_{{0}}a_{{2}}+1/2\, \left( \lambda\,{B_{{1}}}^{2}-\lambda\,{B_{{2}}}^{2}-{\frac {{\mu}^{2}}{\lambda}} \right) {\alpha_{{1}}}^{2}\lambda\,a_{{1}}\\\mbox {{\tt `\#mo("1",mathcolor="white")`}}+5\, \left( -2\, \left( \lambda\,{B_{{1}}}^{2}-\lambda\,{B_{{2}}}^{2}-{\frac {{\mu}^{2}}{\lambda}} \right) \lambda+4\,{\mu}^{2} \right) {\alpha_{{1}}}^{4}\alpha_{{0}}a_{{4}}+10\, \left( \lambda\,{B_{{1}}}^{2}-\lambda\,{B_{{2}}}^{2}-{\frac {{\mu}^{2}}{\lambda}} \right) {\alpha_{{1}}}^{2}{\alpha_{{0}}}^{3}a_{{4}}\\\mbox {{\tt `\#mo("1",mathcolor="white")`}}+6\, \left( \lambda\,{B_{{1}}}^{2}-\lambda\,{B_{{2}}}^{2}-{\frac {{\mu}^{2}}{\lambda}} \right) {\alpha_{{1}}}^{2}{\alpha_{{0}}}^{2}a_{{3}}-6\,\lambda\,{\beta_{{0}}}^{2}{\alpha_{{1}}}^{2}a_{{3}}\\-2\,\lambda\,{\beta_{{0}}}^{2}a_{{5}}\alpha_{{0}}+6\,\mu\,\beta_{{0}}{\alpha_{{1}}}^{2}a_{{2}}+3\,\mu\,\beta_{{0}}a_{{5}}{\alpha_{{0}}}^{2}+3/2\,\mu\,a_{{1}}\alpha_{{0}}\beta_{{0}}\\-36\, \left( \lambda\,{B_{{1}}}^{2}-\lambda\,{B_{{2}}}^{2}-{\frac {{\mu}^{2}}{\lambda}} \right) a_{{1}}{\alpha_{{1}}}^{2}-36\,a_{{1}}{\beta_{{0}}}^{2}=0\\"$

(4)

Download split.mw

A procedure to emulate operator D...

Answered: mmcdara 7072

December 15 2024

1 15

(and nothing more)

restart
alias(F=F(x, t), G=G(x, t)):
                             
with(PDEtools):
undeclare(prime):

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:

An_idea.mw

By the way: I just sent you HERE a reply to your Better simplification and latex export question (LaTeX issue only)

A proposal...

Answered: mmcdara 7072

December 13 2024

2 2

As the b's only intervene through their square I advice you to replace b[i]^2 by B[i] for i=1..3.
Once done I propose you to solve a minimization problem instead of a non linear system.

Let Srw the system after the replacements b[i]^2 = B[i] for i=1..3.
A simple objective function is

J := add(lhs~(Srw)^~2):

Minimize it with suitable consraints.
Here I used four different sets of constrtaints (make your choice or define another set).

>

restart:

>	Digits := 10:

>

L := 1:

>

N := 3:

>	alpha := 1:

>	xexact := t -> t^sqrt(2) + t^sqrt(3):

>	f := simplify(fracdiff(t^sqrt(2), t, alpha)) + simplify(fracdiff(t^sqrt(3), t, alpha)):

>	f := unapply(f, t):

>	xapp := a[0] + sum(a[j]*t^sum(b[i]^2, i = 1 .. j), j = 1 .. N):

>	xapp := unapply(xapp, t):

>	xfrac := sum(a[jj]simplify(GAMMA(sum(b[ii]^2, ii = 1 .. jj) + 1)/GAMMA(sum(b[ii]^2, ii = 1 .. jj) + 1 - alpha))t^(sum(b[ii]^2, ii = 1 .. jj) - alpha), jj = 1 .. N):

>	xfrac := unapply(xfrac, t):

>	xfrac1 := sum(a[jj]simplify(sum(b[ii]^2, ii = 1 .. jj)^(alpha + 1)/(sum(b[ii]^2, ii = 1 .. jj) - alpha))t^(sum(b[ii]^2, ii = 1 .. jj) - alpha), jj = 1 .. N):

>	xfrac1 := unapply(xfrac1, t):

>	S1 := {seq(evalf(xfrac(k/(2N)L)) - evalf(f(k/(2N)L)) = 0, k = 1 .. 2*N)}:

>	S2 := {xapp(0) = 0}:

>	S := S1 union S2:

>

# As the b's only intervene through their square I advice you to replace b[i]^2 by B[i]
# for i=1..3

U   := convert(indets(S, name), list);
bs := select(has, U, b):
Bs := [B[1], B[2], B[3]]:

Srw := eval(S, bs^~2 =~ Bs):   # Set b[i]^2 = B[i]
Urw := eval(U, bs =~ Bs);      # redefine the unknowns

NU := numelems(Urw);

Bp, ap := selectremove(has, Urw, B);

(1)

>	J := add(lhs~(Srw)^~2): # An example of constraints Cstr := { op(Bp >=~ 1), op(ap >=~ 0) }; opt := Optimization:-NLPSolve(J, Cstr); Check := eval(Srw, opt[2]);

[0.659255656578767127e-5, [B[1] = HFloat(1.5425157070011413), B[2] = HFloat(1.0), B[3] = HFloat(1.0826225647977028), a[0] = HFloat(1.865321239496017e-8), a[1] = HFloat(1.9338613155391453), a[2] = HFloat(0.06456529072675758), a[3] = HFloat(0.0)]]

(2)

>	# Another example of constraints Cstr := { B[1] >= 1, B[1]+B[2] >= 1, B[1]+B[2]+B[3] >= 1, op(ap >=~ 0) }; opt := Optimization:-NLPSolve(J, Cstr); Check := eval(Srw, opt[2]);

[0.381760959681121170e-9, [B[1] = HFloat(1.3759235819725324), B[2] = HFloat(0.2077786091095173), B[3] = HFloat(0.18770027396583533), a[0] = HFloat(0.0), a[1] = HFloat(0.6663691594756461), a[2] = HFloat(0.7090748930594262), a[3] = HFloat(0.6246017616247141)]]

(3)

>	# Still another example of constraints Cstr := { B[1] >= 1, B[1]+B[2] >= 1, B[1]+B[2]+B[3] >= 1 }; opt := Optimization:-NLPSolve(J, Cstr); Check := eval(Srw, opt[2]);

[0.548511110584934342e-9, [B[1] = HFloat(1.4451892195800993), B[2] = HFloat(0.34287300100323503), B[3] = HFloat(0.968291120547419), a[0] = HFloat(-2.1900898501097005e-10), a[1] = HFloat(1.2420373086522076), a[2] = HFloat(0.7618141395655179), a[3] = HFloat(-0.003951334337020238)]]

{HFloat(-1.380083329127757e-5) = 0, HFloat(-5.632758017348749e-6) = 0, HFloat(-2.8410628478692246e-6) = 0, HFloat(-2.6872804892441593e-6) = 0, HFloat(-2.1900898501097005e-10) = 0, HFloat(1.194773929169557e-5) = 0, HFloat(1.2972222608542694e-5) = 0}

(4)

>	# And the last constraints example Cstr := { }; opt := Optimization:-NLPSolve(J, Cstr); Check := eval(Srw, opt[2]);

[0.382589102773216340e-7, [B[1] = HFloat(1.618800458856752), B[2] = HFloat(0.08422433693447466), B[3] = HFloat(-0.3287342847332223), a[0] = HFloat(-8.210928420287353e-8), a[1] = HFloat(0.08572187017950603), a[2] = HFloat(1.1448863143053623), a[3] = HFloat(0.7695969717185049)]]

(5)

Download Proposal.mw

Wiener and modulated Wiener processes on...

Answered: mmcdara 7072

December 12 2024

2 0

>

restart:

>	with(Finance):

>	W := WienerProcess():

>	P := PathPlot(W(t), t = 0..3, timesteps = 50, replications = 1): P;

>	W__points := plottools:-getdata(P)[1, -1]

$`#msub(mi("W"),mi("points"))` := Vector(4, {(1) = ` 51 x 2 `*Matrix, (2) = `Data Type: `*float[8], (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})$

(1)

>	f := (t, lambda) -> exp(lambda*t)

(2)

>	Lambda := -1: Z__points := `<\|>`(W__points[..,1], W__points[..,2] *~ f~(W__points[..,1], Lambda))

$`#msub(mi("Z"),mi("points"))` := Vector(4, {(1) = ` 51 x 2 `*Matrix, (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})$

(3)

>	plot([W__points, Z__points], color=[blue, red], legend=[typeset('W'(t)), typeset('W'(t)exp(Lambdat))])

>

Download WienerProcess.mw

Option 1restartgo2 := proc(i) plots:-d...

Answered: mmcdara 7072

December 09 2024

0 0

In case your Maple version does not have Statistics:-HeatMap (which is unnecessary given you have only two colors), I propose you two alternatives based upon plots:-sparsematrixplot

Option 1 (the simplest... IMO)

restart:
go2 := proc() 
  plots:-sparsematrixplot(LinearAlgebra:-RandomMatrix(50, generator = 0 .. 1), axes = none)
end proc:

plots:-animate(go2, [], i = 1 .. 10, frames=10, background=gray, title="")

Option 2 (using Explore the way acer did)

restart
go2 := proc(i) 
  plots:-display(
    plots:-sparsematrixplot(LinearAlgebra:-RandomMatrix(50, generator = 0 .. 1), axes = none)
    , plottools:-rectangle([0.5$2], [50.5$2], color=gray, style=polygon)
  ) 
end proc:

Explore(go2(i), parameters=[i = 1 .. 10], animate, numframes = 10)

By defining a rule?...

Answered: mmcdara 7072

November 19 2024

0 3

MyRule := cos = (theta -> sqrt(1-sin(theta)^2));

Then, e being some expression invoking sines and cosines, remove these latters by executing

eval(e, MyRule)

But you must be more detailed and, by the way, providing an example would help us delivering a better solution.

how to get curve fitting parameters righ...

Answered: mmcdara 7072

November 19 2024

1 4

... by chosing a good model (and this is not a joke !).

PART ONE: OUT OF ANY PHYSICAL CONSIDERATIONS

Looking to your experimental points suggest a fisrst order approximation could be a "capacitor charge curve".
So a model of the form:

(1-exp(-beta*t) / (1-exp(-beta))    (with beta > 0)

Using this extremely simple model already gives a pretty good result provided the addition of a scaling parameter

CapacitorCharge := t ->  scaling * (1 -exp(-beta*t)) / (1-exp(-beta)):

obj := add( ( CapacitorCharge~(c_time) - c_strain )^~2 ):

opt := NLPSolve(obj, {beta >= 1e-6, scaling >= 1e-6 });

 [0.00306881881374611830, [beta = 1.7599456583228286), scaling = 0.16882154828021012]]

display(
   ScatterPlot(c_time, c_strain, symbol=circle, color=blue),
   plot(eval(CapacitorCharge(t), opt[2]), t=0.00..max(c_time), color=black)
)

It's up to you to add more components with the same generic model.
For instance

CapacitorCharge := (t, s, b) ->  s * (1 -exp(-b*t)) / (1-exp(-b)):

TwoComponentModel := unapply(CapacitorCharge(t, s1, b1)  + CapacitorCharge(t, s2, b2), t):

obj := add( ( TwoComponentModel~(c_time) - c_strain )^~2 ):

eps := 1e-6:
opt := Optimization:-NLPSolve(
              obj
              , {
                   s1 >= eps
                   , s2 >= 2*s1  # illustrative example
                   , b1 >= eps
                   , b2 >= 2*b1  # illustrative example
                }
             ):

plots:-display(
   Statistics:-ScatterPlot(c_time, c_strain, symbol=circle, color=blue)
   , plot(eval(TwoComponentModel(t), opt[2]), t=0.00..max(c_time), color=black)
)

Why did you fail?
I suggest you to begin modifying your creep_strain procedure this way

creep_strain := proc(t, n)
  epsilon[0]*(1 + alpha[theta]*add(-(B[i]*(-beta[i]*t + exp(-beta[i] *t) -1))/beta[i],i=1..n))
end proc;

and use Explore to see what happens, in particular how your model (and already its single component version that creep_strain(t, 1) represents) looks almost linear in the range t=0..5

# My Maple 2015 doesn't accept indexed parameters:

One_component_model := eval(creep_strain(t, 1), {B[1] = p, beta[1] = q});

Explore( plot(One_component_model, t=0..5),  parameters=[p=1e-3 .. 1, q=1e-4..1e-2]);

This plot is also quite informative.

plots:-display(
   Statistics:-ScatterPlot(c_time, c_strain, symbol=circle, color=blue)
   , plot(eval(CapacitorCharge(t), opt[2]), t=0.00..max(c_time), color=black)
   , plot(eval(creep_strain(t, 1), eval(opt[2], {scaling=B[1], beta=beta[1]})), t=0.00..max(c_time), color=red)
)

PART TWO: ACCOUNTING FOR PHYSICAL CONSIDERATIONS

Your data come very likely from experiments on some viscoelastic material whose behavior you have inferred could be modeled by your creep_strain function.
As it seems not to be the case, this mean that your hypothesis is not coroborated by the observations, so use another viscoelastic model.
For instance, a simple Kelvin-Voigt model gives this:

(the # The bounds have to be adjusted trick prevents getting identical components, an alternative possibiily using a penalized objective function is given further)

>

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)

>

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

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

obj := add( ( NComponentModel(1)~(c_time) - c_strain )^~2 ):

eps := 1e-6:
opt1 := Optimization:-NLPSolve( obj, {s1 >= eps , b1 >= eps});

p1 := plot(eval(NComponentModel(1)(t), opt1[2]), t=0.00..max(c_time), color=red, legend="1 component"):

(5)

>

obj := add( ( NComponentModel(2)~(c_time) - c_strain )^~2 ):

# The bounds have to be adjusted
opt2 := Optimization:-NLPSolve( obj, {s1 >= eps , b1 >= eps, seq(s||i >= i*s1, i=2..2), seq(b||i <= b1/i, i=2..2)}, iterationlimit=1000);

p2 := plot(eval(NComponentModel(2)(t), opt2[2]), t=0.00..max(c_time), color=gold, legend="2 component"):

(6)

>

obj := add( ( NComponentModel(3)~(c_time) - c_strain )^~2 ):

# The bounds have to be adjusted
opt3 := Optimization:-NLPSolve( obj, {s1 >= eps , b1 >= eps, seq(s||i >= i*s1, i=2..3), seq(b||i <= b1/i, i=2..3)}, iterationlimit=1000);

p3 := plot(eval(NComponentModel(3)(t), opt3[2]), t=0.00..max(c_time), color=green, legend="3 component"):

(7)

>	plots:-display( Statistics:-ScatterPlot(c_time, c_strain, symbol=circle) , p1, p2, p3 )

>

Download Kelvin_Voigt_mmcdara.mw

Here is the penalized version: Kelvin_Voigt_penalized_mmcdara.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

7072 Reputation

18 Badges

MaplePrimes Activity

These are answers submitted by mmcdara

output=Array...

Another way...

For information only...

Use similar triangles...

A less brute force...

Here it is...

'known' option...

Do you mean something like this?...

Nice (?) splitting of a latex expression...

A procedure to emulate operator D...

A proposal...

Wiener and modulated Wiener processes on...

Option 1restartgo2 := proc(i) plots:-d...

By defining a rule?...

how to get curve fitting parameters righ...

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