Advanced Search

tomleslie

13876 Reputation

20 Badges

15 years, 163 days
Contact tomleslie

MaplePrimes Activity


These are answers submitted by tomleslie

Difficult to work out...

tomleslie 13876
May 11 2022
1 0

exactly what your problem is, but may be the "toy" examle in the attached will help?

> 

  restart:
> 

#
# A "toy" example
#
  p1:=plot(Heaviside(x-5)-Heaviside(x-10), x=0..15):
#
# Reflect the above about the line x=0 (the y-axis)
#
  p2:=plottools:-reflect( p1,
                          [[0,0],[0,1]]
                        ):
#
# Display both plots
#
  plots:-display( [p1,p2],
                  color=[red, blue]
                );

 
> 

 


 

Download rplot.mw

Actually no...

tomleslie 13876
May 08 2022
2 3

Because there are several alternative if they are applied correctly - see the attached

> 

  restart;
  ScientificConstants:-GetValue(ScientificConstants:-Constant(g));

9.80665

 (1)
> 

  restart:
  with(ScientificConstants):
  gg:=Constant(g):
  GetValue(gg);
  

9.80665

 (2)
> 

  restart:
  with(ScientificConstants):
  GetValue(Constant(g));

9.80665

 (3)
> 

 

Download SC.mw

 

Maybe...

tomleslie 13876
May 08 2022
2 0

you are after polynomials of a slightly different form -like those in the attached perhaps?

> 

  P:=(x,y,N)->sum(sum(a[i,j]*x^i*y^j,i=0..N), j=0..N);
  P(x,y,3);

proc (x, y, N) options operator, arrow; sum(sum(a[i, j]*x^i*y^j, i = 0 .. N), j = 0 .. N) end proc

  

x^3*y^3*a[3, 3]+x^3*y^2*a[3, 2]+x^2*y^3*a[2, 3]+x^3*y*a[3, 1]+x^2*y^2*a[2, 2]+x*y^3*a[1, 3]+x^3*a[3, 0]+x^2*y*a[2, 1]+x*y^2*a[1, 2]+y^3*a[0, 3]+x^2*a[2, 0]+x*y*a[1, 1]+y^2*a[0, 2]+x*a[1, 0]+y*a[0, 1]+a[0, 0]

 (1)
> 

#
# First derivative wrt x of general polynomial
#
  simplify( diff(P(x,y,N), x) );
#
# First derivative wrt x of polynomial with N=3
#
  simplify( diff(P(x,y,3), x) );
#
# Second derivative wrt x of general polynomial
#
  simplify( diff(P(x,y,N), x$2) );
#
# Second derivative wrt x of polynomial with N=5
#
  simplify( diff(P(x,y,3), x$2) );
#
# First derivative wrt y of general polynomial
#
  simplify( diff(P(x,y,N), y) );
#
# First derivative wrt y of polynomial with N=3
#
  simplify( diff(P(x,y,3), y) );
#
# Second derivative wrt y of general polynomial
#
  simplify( diff(P(x,y,N), y$2) );
#
# Second derivative wrt x of polynomial with N=3
#
  simplify( diff(P(x,y,3), y$2) );

sum(y^j*(sum(a[i, j]*x^(i-1)*i, i = 0 .. N)), j = 0 .. N)

  

(3*x^2*a[3, 3]+2*x*a[2, 3]+a[1, 3])*y^3+(3*x^2*a[3, 2]+2*x*a[2, 2]+a[1, 2])*y^2+(3*x^2*a[3, 1]+2*x*a[2, 1]+a[1, 1])*y+3*x^2*a[3, 0]+2*x*a[2, 0]+a[1, 0]

  

sum(y^j*(sum(a[i, j]*x^(i-2)*i*(i-1), i = 0 .. N)), j = 0 .. N)

  

(6*x*a[3, 3]+2*a[2, 3])*y^3+(6*x*a[3, 2]+2*a[2, 2])*y^2+(6*x*a[3, 1]+2*a[2, 1])*y+6*x*a[3, 0]+2*a[2, 0]

  

sum(y^(j-1)*j*(sum(a[i, j]*x^i, i = 0 .. N)), j = 0 .. N)

  

(3*y^2*a[3, 3]+2*y*a[3, 2]+a[3, 1])*x^3+(3*y^2*a[2, 3]+2*y*a[2, 2]+a[2, 1])*x^2+(3*y^2*a[1, 3]+2*y*a[1, 2]+a[1, 1])*x+3*y^2*a[0, 3]+2*y*a[0, 2]+a[0, 1]

  

sum(j*(j-1)*y^(j-2)*(sum(a[i, j]*x^i, i = 0 .. N)), j = 0 .. N)

  

(6*y*a[3, 3]+2*a[3, 2])*x^3+(6*y*a[2, 3]+2*a[2, 2])*x^2+(6*y*a[1, 3]+2*a[1, 2])*x+6*y*a[0, 3]+2*a[0, 2]

 (2)
> 

 

Download polydiff2.mw

If I interpret...

tomleslie 13876
May 08 2022
3 4

your question correctly, then the attached ought to contain what you want.

> 

  restart;
  P:=(x,y,N)->sum( alpha[i]*x^i*y^(N-i), i=0..N);
  P(x,y,5);
#
# First derivative wrt x of general polynomial
#
  simplify( diff(P(x,y,N), x) );
#
# First derivative wrt x of polynomial with N=5
#
  simplify( diff(P(x,y,5), x) );
#
# Second derivative wrt x of general polynomial
#
  simplify( diff(P(x,y,N), x$2) );
#
# Second derivative wrt x of polynomial with N=5
#
  simplify( diff(P(x,y,5), x$2) );
#
# First derivative wrt y of general polynomial
#
  simplify( diff(P(x,y,N), y) );
#
# First derivative wrt y of polynomial with N=5
#
  simplify( diff(P(x,y,5), y) );
#
# Second derivative wrt y of general polynomial
#
  simplify( diff(P(x,y,N), y$2) );
#
# Second derivative wrt x of polynomial with N=5
#
  simplify( diff(P(x,y,5), y$2) );
  

proc (x, y, N) options operator, arrow; sum(alpha[i]*x^i*y^(N-i), i = 0 .. N) end proc

  

x^5*alpha[5]+x^4*y*alpha[4]+x^3*y^2*alpha[3]+x^2*y^3*alpha[2]+x*y^4*alpha[1]+y^5*alpha[0]

  

sum(alpha[i]*x^(i-1)*i*y^(N-i), i = 0 .. N)

  

5*x^4*alpha[5]+4*x^3*y*alpha[4]+3*x^2*y^2*alpha[3]+2*x*y^3*alpha[2]+y^4*alpha[1]

  

sum(alpha[i]*x^(i-2)*i*y^(N-i)*(i-1), i = 0 .. N)

  

20*x^3*alpha[5]+12*x^2*y*alpha[4]+6*x*y^2*alpha[3]+2*y^3*alpha[2]

  

sum(alpha[i]*x^i*y^(N-i-1)*(N-i), i = 0 .. N)

  

x^4*alpha[4]+2*x^3*y*alpha[3]+3*x^2*y^2*alpha[2]+4*x*y^3*alpha[1]+5*y^4*alpha[0]

  

sum(alpha[i]*x^i*y^(N-i-2)*(N-i)*(N-i-1), i = 0 .. N)

  

2*x^3*alpha[3]+6*x^2*y*alpha[2]+12*x*y^2*alpha[1]+20*y^3*alpha[0]

 (1)
> 

 


Download polydiff.mw

Try...

tomleslie 13876
May 07 2022
1 1

the code in the attached

> 

  restart;
> 

  S:= [ seq
        ( add
          ( parse~
            ( StringTools:-Explode
                           ( convert
                             ( convert
                               ( j,
                                 binary
                               ),
                               string
                             )
                           )
            )
          ),
          j=0..10
       )
     ];

[0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2]

 (1)
> 

 


 

Download numOnes.mw

 

The two methods...

tomleslie 13876
May 07 2022
2 1

give very similar answers, but they are different. A few observations

  1. Both methods are "good" ways to get the numeric solution of an ODE system.
  2. Since both are going to get close to the"actual" numeric solution, then they must be close to each other
  3. The ODE system you are using has smoothly varying solutions, which makes things "easy" for a numeric solver. Even methods less sophisticated than those you have tried, would probably get pretty close to the actual solution for this problem

In the attached I have added an execution group which shows the discrepancy between the two methods. Since the calculation in the worksheet is performed with the default setting of Digits (ie 10), one might expect discrepancies due to straightforward floating-point rounding to be ~O(10-9). The discrepancies between your two methods over most of the points are in the range 10-5 to 10-7 That isn't floating point rounding!


 

> 

restart; with(PDEtools, Solve); with(LinearAlgebra); with(plots); Digits := 15; N := 8; t0, tf := 0, 1; Ntt := 10; h := evalf((tf-t0)/(Ntt-1))

0, 1

  

.111111111111111

 (1)
> 

sys := [diff(A[2](t), t) = (4/45)*(1323*Pi^5*A[7](t)+945*Pi^4*A[6](t)+630*Pi^3*A[5](t)+513*Pi^2*A[4](t)+189*Pi*A[3](t)-3626*A[2](t))/Pi^2, diff(A[3](t), t) = -(4/45)*(9744*Pi^5*A[7](t)+6960*Pi^4*A[6](t)+4415*Pi^3*A[5](t)+2784*Pi^2*A[4](t)+1392*Pi*A[3](t)-14544*A[2](t))/Pi^3, diff(A[4](t), t) = (2/9)*(14112*Pi^5*A[7](t)+10215*Pi^4*A[6](t)+6720*Pi^3*A[5](t)+4032*Pi^2*A[4](t)+2016*Pi*A[3](t)-12992*A[2](t))/Pi^4, diff(A[5](t), t) = -(2/45)*(133455*Pi^5*A[7](t)+96000*Pi^4*A[6](t)+64000*Pi^3*A[5](t)+38400*Pi^2*A[4](t)+19200*Pi*A[3](t)-81408*A[2](t))/Pi^5, diff(A[6](t), t) = (4096/45)*(63*Pi^5*A[7](t)+45*Pi^4*A[6](t)+30*Pi^3*A[5](t)+18*Pi^2*A[4](t)+9*Pi*A[3](t)-26*A[2](t))/Pi^6, diff(A[7](t), t) = -(32768/315)*(21*Pi^5*A[7](t)+15*Pi^4*A[6](t)+10*Pi^3*A[5](t)+6*Pi^2*A[4](t)+3*Pi*A[3](t)-6*A[2](t))/Pi^7]; nops(sys)
> 

unknownAs := {seq(A[i](t), i = 2 .. N-1)}
> 

f, diffs := GenerateMatrix(`~`[`-`](`~`[rhs](sys), `~`[lhs](sys)), [op(1 .. N-2, unknownAs)])

Matrix(%id = 36893488148077213564), Vector[column](%id = 36893488148077213444)

 (2)
> 

NULL
> 

Digits := 20
> 

ICS := [A[2](0) = -0.121093412283174e-2, A[3](0) = -.163637479659224, A[4](0) = -0.383881333426826e-2, A[5](0) = 0.110064697379063e-1, A[6](0) = -0.100570735132924e-2, A[7](0) = -0.277528316343890e-4]
> 

``
> 

> 

npts := Ntt; resMat4 := Matrix(N-2, npts); resMat4[1 .. N-2, 1] := Vector(`~`[rhs](ICS)); for j to Ntt-1 do k41 := f.resMat4[() .. (), j]; k42 := f.(resMat4[() .. (), j]+(1/2)*h*k41); k43 := f.(resMat4[() .. (), j]+(1/2)*h*k42); k44 := f.(h*k43+resMat4[() .. (), j]); resMat4[() .. (), j+1] := resMat4[() .. (), j]+h*((1/6)*k41+(1/3)*k42+(1/3)*k43+(1/6)*k44) end do
> 

npts := Ntt; resMatkf := Matrix(N-2, npts); resMatkf[1 .. N-2, 1] := Vector(`~`[rhs](ICS)); for j to Ntt-1 do kf1 := f.resMatkf[() .. (), j]; kf2 := f.(resMatkf[() .. (), j]+(1/4)*h*kf1); kf3 := f.(resMatkf[() .. (), j]+(3/32)*h*kf1+(9/32)*h*kf2); kf4 := f.(resMatkf[() .. (), j]+(1932/2197)*h*kf1-(7200/2197)*h*kf2+(7296/2197)*h*kf3); kf5 := f.(resMatkf[() .. (), j]+(439/216)*h*kf1-8*h*kf2+(3680/513)*h*kf3-(845/4104)*h*kf4); kf6 := f.(resMatkf[() .. (), j]-(8/27)*h*kf1+2*h*kf2-(3544/2565)*h*kf3+(1859/4104)*h*kf4-(11/40)*h*kf5); resMatkf[() .. (), j+1] := resMatkf[() .. (), j]+h*((16/135)*kf1+(6656/12825)*kf3+(28561/56430)*kf4-(9/50)*kf5+(2/55)*kf6) end do
> 

interface(displayprecision = 5); resMat4; resMatkf; resMat4-resMatkf

Matrix(6, 10, {(1, 1) = -0.121e-2, (1, 2) = -0.127e-2, (1, 3) = -0.125e-2, (1, 4) = -0.117e-2, (1, 5) = -0.107e-2, (1, 6) = -0.96e-3, (1, 7) = -0.86e-3, (1, 8) = -0.76e-3, (1, 9) = -0.68e-3, (1, 10) = -0.61e-3, (2, 1) = -.16364, (2, 2) = -.14647, (2, 3) = -.13089, (2, 4) = -.11702, (2, 5) = -.10468, (2, 6) = -0.9366e-1, (2, 7) = -0.8381e-1, (2, 8) = -0.7500e-1, (2, 9) = -0.6712e-1, (2, 10) = -0.6006e-1, (3, 1) = -0.384e-2, (3, 2) = -0.311e-2, (3, 3) = -0.294e-2, (3, 4) = -0.273e-2, (3, 5) = -0.248e-2, (3, 6) = -0.223e-2, (3, 7) = -0.199e-2, (3, 8) = -0.177e-2, (3, 9) = -0.158e-2, (3, 10) = -0.141e-2, (4, 1) = 0.1101e-1, (4, 2) = 0.957e-2, (4, 3) = 0.863e-2, (4, 4) = 0.778e-2, (4, 5) = 0.699e-2, (4, 6) = 0.626e-2, (4, 7) = 0.560e-2, (4, 8) = 0.500e-2, (4, 9) = 0.447e-2, (4, 10) = 0.400e-2, (5, 1) = -0.101e-2, (5, 2) = -0.80e-3, (5, 3) = -0.73e-3, (5, 4) = -0.67e-3, (5, 5) = -0.61e-3, (5, 6) = -0.54e-3, (5, 7) = -0.49e-3, (5, 8) = -0.44e-3, (5, 9) = -0.39e-3, (5, 10) = -0.35e-3, (6, 1) = -0.3e-4, (6, 2) = -0.4e-4, (6, 3) = -0.3e-4, (6, 4) = -0.3e-4, (6, 5) = -0.2e-4, (6, 6) = -0.2e-4, (6, 7) = -0.2e-4, (6, 8) = -0.2e-4, (6, 9) = -0.2e-4, (6, 10) = -0.1e-4})

  

Matrix(6, 10, {(1, 1) = -0.121e-2, (1, 2) = -0.144e-2, (1, 3) = -0.130e-2, (1, 4) = -0.118e-2, (1, 5) = -0.107e-2, (1, 6) = -0.96e-3, (1, 7) = -0.86e-3, (1, 8) = -0.76e-3, (1, 9) = -0.68e-3, (1, 10) = -0.61e-3, (2, 1) = -.16364, (2, 2) = -.14613, (2, 3) = -.13080, (2, 4) = -.11701, (2, 5) = -.10468, (2, 6) = -0.9366e-1, (2, 7) = -0.8381e-1, (2, 8) = -0.7500e-1, (2, 9) = -0.6712e-1, (2, 10) = -0.6006e-1, (3, 1) = -0.384e-2, (3, 2) = -0.343e-2, (3, 3) = -0.302e-2, (3, 4) = -0.274e-2, (3, 5) = -0.248e-2, (3, 6) = -0.223e-2, (3, 7) = -0.199e-2, (3, 8) = -0.177e-2, (3, 9) = -0.158e-2, (3, 10) = -0.141e-2, (4, 1) = 0.1101e-1, (4, 2) = 0.972e-2, (4, 3) = 0.866e-2, (4, 4) = 0.778e-2, (4, 5) = 0.699e-2, (4, 6) = 0.626e-2, (4, 7) = 0.560e-2, (4, 8) = 0.500e-2, (4, 9) = 0.447e-2, (4, 10) = 0.400e-2, (5, 1) = -0.101e-2, (5, 2) = -0.84e-3, (5, 3) = -0.74e-3, (5, 4) = -0.67e-3, (5, 5) = -0.61e-3, (5, 6) = -0.54e-3, (5, 7) = -0.49e-3, (5, 8) = -0.44e-3, (5, 9) = -0.39e-3, (5, 10) = -0.35e-3, (6, 1) = -0.3e-4, (6, 2) = -0.3e-4, (6, 3) = -0.3e-4, (6, 4) = -0.3e-4, (6, 5) = -0.2e-4, (6, 6) = -0.2e-4, (6, 7) = -0.2e-4, (6, 8) = -0.2e-4, (6, 9) = -0.2e-4, (6, 10) = -0.1e-4})

  

Matrix(%id = 36893488148273947940)

 (3)
> 

``


 

Download matDiff.mw

I wouldn't bother...

tomleslie 13876
May 06 2022
2 0

Your code contains multiple errors, which prevent it running successfully in Maple.

Why transfer to Matlab? Maybe you think that transferring  code with multiple errors will suddenly become "functional"??? That's an impossible dream!. Transferring perfectly functional code can be problematic: transferring code with known errors - well you will transfer all the errors (and maybe generate a few more).

I suggest that you read the comments in the final execution group of the attached very carefully. To summarize, you have two main issues

  1. The NLPSolve() uses four constraints, two of which are not valid constraints, because they are simple numbers, not equalities or inequalities: the other two contain multiple variables x1[1], x1[3], x2[1], x2[2], x2[3], which are undefined and do not occur in the expression to be minimized
  2. The final plot() command complains because the expression x1t contains references to the variables x1[1], x1[2], x1[3], none of which have been defined
> 

restart; t1 := time(); with(LinearAlgebra); J := readstat("Please enter integer number J: "); N1 := proc (x) options operator, arrow; piecewise(0 <= x and x <= 1, 1) end proc; N2 := proc (x) options operator, arrow; piecewise(0 <= x and x <= 1, x, 1 < x and x <= 2, 2-x) end proc; N := proc (J, k) options operator, arrow; unapply(N2(2^J*x-k), x) end proc; Phi := proc (J, k) options operator, arrow; evalf((N(J, k))(x))*N1(x) end proc; PhiJ := Vector[column](2^J+1); for k from -1 to 2^J-1 do PhiJ[k+2] := Phi(J, k) end do; P := Matrix(2^J+1, 2^J+1); Map2[proc (i, j) options operator, arrow; evalb(i-j = 1) end proc](proc (x, a) options operator, arrow; x end proc, 1/6, P, inplace); Map2[proc (i, j) options operator, arrow; evalb(j-i = 1) end proc](proc (x, a) options operator, arrow; x end proc, 1/6, P, inplace); Map2[proc (i, j) options operator, arrow; evalb(i = j) end proc](proc (x, a) options operator, arrow; x end proc, 2/3, P, inplace); P[1, 1] := 1/3; P[2^J+1, 2^J+1] := 1/3; P := 2^(-J)*P; E := Matrix(2^J+1, 2^J+1); Map2[proc (i, j) options operator, arrow; evalb(i-j = 1) end proc](proc (x, a) options operator, arrow; x end proc, 1/2, E, inplace); Map2[proc (i, j) options operator, arrow; evalb(j-i = 1) end proc](proc (x, a) options operator, arrow; x end proc, -1/2, E, inplace); E[1, 1] := -1/2; E[2^J+1, 2^J+1] := 1/2; DPhi := E.(1/P); X1 := Vector[column](2^J+1, symbol = x1); X2 := Vector[column](2^J+1, symbol = x2); U := Vector[column](2^J+1, symbol = u); JJ := (1/2)*U^%T.P.U; x1t := X1^%T.PhiJ; x2t := X2^%T.PhiJ; ut := U^%T.PhiJ; for i from 0 to 2^J do PhiJxJ[i+1] := apply(unapply(PhiJ, x), i/2^J) end do; for i to 2^J+1 do eq1[i] := (X1^%T.DPhi-X2^%T).PhiJxJ[i] = 0; eq2[i] := (X2^%T.DPhi-U^%T).PhiJxJ[i] = 0 end do; for i to 2^J+1 do eq3[i] := X1^%T.PhiJxJ[i]-.1, 0 end do; eq1[0] := eval(x1t, x = 0) = 0; eq2[0] := eval(x2t, x = 0)-1 = 0; eq1[2^J+2] := eval(x1t, x = 1) = 0; eq2[2^J+2] := eval(x2t, x = 1) = -1; eqq1 := {seq(eq1[i], i = 0*.2^J+2)}; eqq2 := {seq(eq2[i], i = 0.2^J+2)}; eqq3 := {seq(eq3[i], i = 1.2^J+1)}; eq := `union`(`union`(eqq1, eqq2), eqq3); with(Optimization); S := NLPSolve(JJ, eq); assign(S[2]); uexact := piecewise(0 <= x and x <= .3, (200/9)*x-20/3, .3 <= x and x <= .7, 0, .7 <= x and x <= 1, -(200/9)*x+140/9); x2exact := piecewise(0 <= x and x <= .3, (100/9)*x^2-(20/3)*x+1, .3 <= x and x <= .7, 0, .7 <= x and x <= 1, -(100/9)*x^2+(140/9)*x-49/9); x1exact := piecewise(0 <= x and x <= .3, (100/27)*x^3-(10/3)*x^2+x, .3 <= x and x <= .7, 1/10, .7 <= x and x <= 1, -(100/27)*x^3+(70/9)*x^2-(49/9)*x+37/27); plot([x1exact, x1t], x = 0 .. 1, style = [line, point], legend = ["Exact", "Approximate"], axis = [gridlines = [colour = green, majorlines = 2]], labels = ["t", x[1](t)], labeldirections = ["horizontal", "vertical"])

Error, (in Optimization:-NLPSolve) constraints must be specified as a set or list of equalities and inequalities

  

Warning, expecting only range variable x in expression x1[1]*piecewise(0. <= 2.*x+1. and 2.*x <= 0.,2.*x+1.,0. < 2.*x and 2.*x <= 1.,1.-2.*x)*piecewise(0 <= x and x <= 1,1)+x1[2]*piecewise(0. <= 2.*x and 2.*x <= 1.,2.*x,1. < 2.*x and 2.*x <= 2.,2.-2.*x)*piecewise(0 <= x and x <= 1,1)+x1[3]*piecewise(0. <= 2.*x-1. and 2.*x <= 2.,2.*x-1.,0. < 2.*x-2. and 2.*x <= 3.,3.-2.*x)*piecewise(0 <= x and x <= 1,1) to be plotted but found names [x1[1], x1[2], x1[3]]

  

 
> 

JJ; indets(JJ, name); eq; `union`(`~`[indets](eq, name)[]); x1exact; x1t; indets(x1exact, name); indets(x1t, name)

((1/12)*u[1]+(1/24)*u[2])*u[1]+((1/24)*u[1]+(1/6)*u[2]+(1/24)*u[3])*u[2]+((1/24)*u[2]+(1/12)*u[3])*u[3]

  

{u[1], u[2], u[3]}

  

{0, 1.*x1[3]-.1, -1.000000000*x1[1]+1.000000000*x1[3]-1.000000000*x2[2] = 0, -1.000000000*x2[1]+1.000000000*x2[3]-1.000000000*u[2] = 0}

  

{u[2], x1[1], x1[3], x2[1], x2[2], x2[3]}

  

piecewise(0 <= x and x <= .3, (100/27)*x^3-(10/3)*x^2+x, .3 <= x and x <= .7, 1/10, .7 <= x and x <= 1, -(100/27)*x^3+(70/9)*x^2-(49/9)*x+37/27)

  

x1[1]*piecewise(0. <= 2.*x+1. and 2.*x <= 0., 2.*x+1., 0. < 2.*x and 2.*x <= 1., 1.-2.*x)*piecewise(0 <= x and x <= 1, 1)+x1[2]*piecewise(0. <= 2.*x and 2.*x <= 1., 2.*x, 1. < 2.*x and 2.*x <= 2., 2.-2.*x)*piecewise(0 <= x and x <= 1, 1)+x1[3]*piecewise(0. <= 2.*x-1. and 2.*x <= 2., 2.*x-1., 0. < 2.*x-2. and 2.*x <= 3., 3.-2.*x)*piecewise(0 <= x and x <= 1, 1)

  

{x}

  

{x, x1[1], x1[2], x1[3]}

 (1)

 

Download badCode.mw

Use surd()...

tomleslie 13876
May 06 2022
2 1

as in the attached (because when x is negative, you don't want the principal root -ti's complex)

> 

  restart;
  plot(surd(x,3), x=-2..2);

 
> 

 

Download surd.mw

 

I'm (probably?) misunderstanding this qu...

tomleslie 13876
May 03 2022
2 0

but if you want a series solution to a differentail equation in Maple, then you call dsolve() with the 'series' option, as in the attached.

If this is not what you want then I think you are going to have to explain your requirement more clearly

> 

  restart;
> 

  ode := diff(u(t), t) = f(t, u(t));
#
# Second order series solution (NB By default
# Order =6
#
  Order := 2:
  dsolve(ode, u(t), series);
#
# Fourth order series solution
#
  Order := 4:
  dsolve(ode, u(t), series);

diff(u(t), t) = f(t, u(t))

  

u(t) = series(u(0)+f(0, u(0))*t+O(t^2),t,2)

  

u(t) = series(u(0)+f(0, u(0))*t+((1/2)*(D[1](f))(0, u(0))+(1/2)*(D[2](f))(0, u(0))*f(0, u(0)))*t^2+((1/6)*(D[2](f))(0, u(0))^2*f(0, u(0))+(1/6)*f(0, u(0))^2*(D[2, 2](f))(0, u(0))+(1/6)*(D[2](f))(0, u(0))*(D[1](f))(0, u(0))+(1/3)*(D[1, 2](f))(0, u(0))*f(0, u(0))+(1/6)*(D[1, 1](f))(0, u(0)))*t^3+O(t^4),t,4)

 (1)
> 

``

NULL

Download odeSer.mw

Or recursively?...

tomleslie 13876
May 03 2022
2 0

as in the attached

> 

  restart:
  bj:= proc( func, q, L1, ind:=1)
             if   ind < numelems(L1)
             then return bj( func,
                             func( q, L1[ind] )-v[ind],
                             L1,
                             ind+1
                           );
                  else return func( q,L1[ind] )-v[ind];
                  fi;
       end proc:
  v:= [p,q,r,s];
  bj( F, y, [a,b,c,d]);

[p, q, r, s]

  

F(F(F(F(y, a)-p, b)-q, c)-r, d)-s

 (1)
> 

  v:='v';
  bj( F, y, [a,b,c,d]);
  eval(%, v=[p, q, r, s]);

v

  

F(F(F(F(y, a)-v[1], b)-v[2], c)-v[3], d)-v[4]

  

F(F(F(F(y, a)-p, b)-q, c)-r, d)-s

 (2)
> 

 


Download oddProc.mw

Hmm...

tomleslie 13876
May 02 2022
3 2

In future please supply a worksheet using the big green up-arrow in the Mapleprimes toolbar. Few responders here are prepared to spend the time typing in your worksheet in order to investigate your problem.

Two points

  1. I cannot reproduce the first problem which you report about "The above gives internal error." In the attached, I have tried you example with and without an expansion point, with and without boundary conditions, with and without uneval quotes around the word 'series'. Since you haven't provided a worksheet showing this "internal error", I'm going to pass on this (for now).
  2. Your second example of what happens when the series solution is complex is more interesting. I think I can see why the odetest() command *might be* failing. Conside the scenario
    1. For some ODE, one obtains a series solution up to (say) O(x6) . NB assume no complex powers of the independent variable
    2. Now this is not an "exact" solution so odetest() presumably decides that it is a "valid" series solution, provided that (on back-subsitution into the ODE), any discrepancies in the ode or initiial conditions are O(x7 or higher)
    3. Suppose now you have a new problem, whose series solution does contain complex powers of the independent variable.
    4. In  the attached I have back-substituted the series solution for your second example into the associated ODE "manually". You can see from the output of the final execution group that such back-substitution contains terms such as x^6*x^(1/2 - sqrt(3)*I/2) - ie the independent variable has a complex exponent. How does one then determine whether any discrepancies in the ode or initiial conditions are O(x7 or higher) - the question doesn't make sense (to me) at all!
> 

  restart;
  interface(version);

`Standard Worksheet Interface, Maple 2022.0, Windows 7, March 8 2022 Build ID 1599809`

 (1)
> 

#
# Example from Maple help just to check everything is working
# (more or less)
#
  ODE := [diff(y(x),x,x)+diff(y(x),x)^2=0, y(a)=0, D(y)(a)=1];
  sol := dsolve( ODE, y(x), type='series');
  odetest(sol, ODE, series);

[diff(diff(y(x), x), x)+(diff(y(x), x))^2 = 0, y(a) = 0, (D(y))(a) = 1]

  

y(x) = series(x-a-(1/2)*(x-a)^2+(1/3)*(x-a)^3-(1/4)*(x-a)^4+(1/5)*(x-a)^5+O((x-a)^6),x = a,6)

  

[0, 0, 0]

 (2)
> 

#
# OP's first example (with no initial conditions and no
# expansion point)
#
# Works as expected
#
  restart;
  ODE:= diff(y(x), x$2)=(3*x^2+c)*diff(y(x),x)+((3-b)*x-a)*y(x);
  Order:= 5;
  sol := dsolve( ODE, y(x), type=series);
  odetest(sol, ODE, series);

diff(diff(y(x), x), x) = (3*x^2+c)*(diff(y(x), x))+((3-b)*x-a)*y(x)

  

5

  

y(x) = series(y(0)+(D(y))(0)*x+((1/2)*c*(D(y))(0)-(1/2)*a*y(0))*x^2+((1/6)*(D(y))(0)*c^2-(1/6)*y(0)*a*c-(1/6)*a*(D(y))(0)-(1/6)*y(0)*b+(1/2)*y(0))*x^3+((1/24)*(D(y))(0)*c^3-(1/24)*y(0)*a*c^2-(1/12)*c*a*(D(y))(0)-(1/24)*c*y(0)*b+(1/8)*c*y(0)+(1/2)*(D(y))(0)+(1/24)*a^2*y(0)-(1/12)*(D(y))(0)*b)*x^4+O(x^5),x,5)

  

0

 (3)
> 

#
# OP's first example (with no initial conditions, but with
# an expansion point)
#
# Works as expected
#
  restart;
  ODE:= diff(y(x), x$2)=(3*x^2+c)*diff(y(x),x)+((3-b)*x-a)*y(x);
  Order:= 5;
  sol := dsolve( ODE, y(x), type=series);
  odetest(sol, ODE, series, point=0);

diff(diff(y(x), x), x) = (3*x^2+c)*(diff(y(x), x))+((3-b)*x-a)*y(x)

  

5

  

y(x) = series(y(0)+(D(y))(0)*x+((1/2)*c*(D(y))(0)-(1/2)*a*y(0))*x^2+((1/6)*(D(y))(0)*c^2-(1/6)*y(0)*a*c-(1/6)*a*(D(y))(0)-(1/6)*y(0)*b+(1/2)*y(0))*x^3+((1/24)*(D(y))(0)*c^3-(1/24)*y(0)*a*c^2-(1/12)*c*a*(D(y))(0)-(1/24)*c*y(0)*b+(1/8)*c*y(0)+(1/2)*(D(y))(0)+(1/24)*a^2*y(0)-(1/12)*(D(y))(0)*b)*x^4+O(x^5),x,5)

  

0

 (4)
> 

#
# OP's first example (with initial conditions)
#
# Works as expected
#
  restart;
  ODE:= [ diff(y(x), x$2)=(3*x^2+c)*diff(y(x),x)+((3-b)*x-a)*y(x), y(0)=1,  D(y)(0)=0];
  Order:= 5;
  sol := dsolve( ODE, y(x), type=series);
  odetest(sol, ODE, series);

[diff(diff(y(x), x), x) = (3*x^2+c)*(diff(y(x), x))+((3-b)*x-a)*y(x), y(0) = 1, (D(y))(0) = 0]

  

5

  

y(x) = series(1-((1/2)*a)*x^2+(-(1/6)*c*a-(1/6)*b+1/2)*x^3+(-(1/24)*c^2*a-(1/24)*c*b+(1/8)*c+(1/24)*a^2)*x^4+O(x^5),x,5)

  

[0, 0, 0]

 (5)
> 

#
# OP's second example - series has complex coefficients
#
  restart;
  ODE:=  x^2*diff(y(x), x$2)+x^2*diff(y(x),x)+y(x)=0;
#
# First check exact solution (assuming one can be obtained
#
  sol:= dsolve(ODE);
  odetest(sol, ODE);
#
# Well that worked!
#
# Now look at series solution
#
  Order:= 6;
  sol := simplify(dsolve( ODE, y(x), type=series));
  odetest(sol, ODE, series, point=0);
#
# And that didn't. Try inserting the series solution
# into the ODE "manually"

x^2*(diff(diff(y(x), x), x))+x^2*(diff(y(x), x))+y(x) = 0

  

y(x) = _C1*exp(-(1/2)*x)*x^(1/2)*BesselI(((1/2)*I)*3^(1/2), (1/2)*x)+_C2*exp(-(1/2)*x)*x^(1/2)*BesselK(((1/2)*I)*3^(1/2), (1/2)*x)

  

0

  

6

  

y(x) = _C1*x^(1/2-((1/2)*I)*3^(1/2))*(series(1-(1/2)*x+((I*3^(1/2)-3)/((8*I)*3^(1/2)-16))*x^2+((-I*3^(1/2)+5)/((48*I)*3^(1/2)-96))*x^3+((1/384)*(I*3^(1/2)-5)*(I*3^(1/2)-7)/((I*3^(1/2)-4)*(I*3^(1/2)-2)))*x^4-((1/3840)*(I*3^(1/2)-7)*(I*3^(1/2)-9)/((I*3^(1/2)-4)*(I*3^(1/2)-2)))*x^5+O(x^6),x,6))+_C2*x^(1/2+((1/2)*I)*3^(1/2))*(series(1-(1/2)*x+((I*3^(1/2)+3)/((8*I)*3^(1/2)+16))*x^2+((-I*3^(1/2)-5)/((48*I)*3^(1/2)+96))*x^3+((1/384)*(I*3^(1/2)+5)*(I*3^(1/2)+7)/((I*3^(1/2)+4)*(I*3^(1/2)+2)))*x^4-((1/3840)*(I*3^(1/2)+7)*(I*3^(1/2)+9)/((I*3^(1/2)+4)*(I*3^(1/2)+2)))*x^5+O(x^6),x,6))

  

Error, (in odetest/series) complex argument to max/min: 13/2-1/2*I*3^(1/2)

  
> 

simplify(eval(ODE, convert(sol, polynom)));

-(89/36480)*(_C1*(I*3^(1/2)+261/89)*x^(1/2-((1/2)*I)*3^(1/2))-_C2*x^(1/2+((1/2)*I)*3^(1/2))*(I*3^(1/2)-261/89))*x^6 = 0

 (6)
> 

 

> 

 

Download testODE.mw

Why...

tomleslie 13876
May 02 2022
1 1

use the Jacobi method to solve a simple matrix equation?

Why not just use the in-built LinearSolve() command?

In the attached worksheet I produce solution both ways - with the LinearSolve() command and with "Jacobi iteration". As expected, hese produce the same answer (within the specified tolerance), but neither agrees with the exact solution. The "shape" of the curves produced is spot-on, but the "scale" is incorrect.

I suspect that the discrepancy between the numerical and exact solutions is due to the fact that the matrix 'A' has to be scaled by (some power of?) the step size being used. I'm a bit pushed for time right now, so I'll leave this last issue for you to sort out!

> 

  restart;
  with(LinearAlgebra):
  ode := diff(u(x), x, x) = -2*x^2 + 1;
  Bcs := u(0) = 0, u(1) = 0;
  u_exact := pdsolve({Bcs, ode});
  p1:= plot( rhs(u_exact), x = 0 .. 1, color = red);

diff(diff(u(x), x), x) = -2*x^2+1

  

u(0) = 0, u(1) = 0

  

u(x) = -(1/6)*x*(x^3-3*x+2)

  

 
> 

#
# Set up some system parameters
#
  tol := 0.0001:
  llim := 0:
  ulim := 1:
  f := x -> 1 - 2*x^2:
  n := 128:
#
# Set up required matrices
#
  A:= Matrix
      ( n-1,
        n-1,
        (i, j) -> `if`( i = j,
                        -2,
                        `if`( abs(i - j) = 1,
                              1,
                              0
                            )
                      )
             ):
  b:= < seq
        ( evalf( f(j/n) ),
          j = 1 .. n-1
        )
      >:
> 

#
# Calculation using Maple's built-in LinearSolve()
#
  ls:= LinearSolve(A, b):
  p2:= plot
       ( [ [ llim, 0 ],
           seq
           ( [ j/n, ls[j] ],
             j = 1 .. n-1
           ),
           [ ulim, 0 ]
         ],
         color=blue,
         style=point
       );

 
> 

#
# Same calculation using "Jacobi iteration" (takes longer
# but appears to give the same answer as LinearSolve()
# above
#
# First set up required matrices
#
# Initial guess at result
#
  Xold:= Vector(n-1, 0.0):
#
# Inverse of the matrix of diagonal elements of A
#
   Dinv:= MatrixInverse
          ( Matrix
            ( op(1,A),
              (i,j)-> `if`
                      ( i=j,
                        A[j,j],
                        0
                      )
            )
          ):
#
# Matrix of off-diagonal elements of A
#
  R:= Matrix
      ( op(1,A),
        (i,j)-> `if`
                ( abs(i-j)=1,
                  A[i,j],
                  0
                )
      ):
#
# Now run the iterative calculation - this is the slow bit!
#
  relerr:=1:
  while relerr > tol do
        Xnew:= Dinv.(b-R.Xold);
        relerr:= norm(Xnew-Xold);
        Xold:= Xnew;
  od:
  p3:= plot
       ( [ [ llim, 0 ],
           seq
           ( [j/n, Xnew[j]],
             j = 1 .. n-1
           ),
           [ ulim, 0 ]
         ],
         color=blue,
         style=point
       );

 
> 

 

> 

NULL

NULL

Download jac.mw

So...

tomleslie 13876
May 01 2022
1 0

if I guess your ineptly coded expression correctly, I would get something like the attached?

> 

  restart;
  expr:=exp(x*cos(x));
#
# OP wants to know the area shaded in the following
#
  plots:-shadebetween( 0, expr, x=-Pi/2..Pi/2);
#
# Which would be
#
  evalf(Int( expr, x=-Pi/2..Pi/2));
#
# Idle curiosity - does the above have some kind of
# analytic expression?
#
  identify(%);

exp(x*cos(x))

  

 

3.399914257

  

-(5/9)*3^(1/2)+(8/7)*exp(1)+(8/7)*ln(3)

 (1)
> 

 

Download areaReg.mw

lots of typos...

tomleslie 13876
April 29 2022
2 1

mainly because of your use of whitespace in 2D-input mode. Also some logical errors eg u__g defined as an expression by used as a function?: six format codes in a printf() statement but only four variabless to print.

So after a fair bit of guessing on my part about what you intended, I came up with the attached.

NULL

> 

  restart;
  n := 3;
  X := [seq(k/n, k = 0 .. n)];
  u_g := x->x^2 + x + 1;

3

  

[0, 1/3, 2/3, 1]

  

proc (x) options operator, arrow; x^2+x+1 end proc

 (1)
> 

  printf("        x      u_exact     u_g      error \n");
  for k to numelems(X) do
      err[k] := abs(u_g(X[k]) - exp(X[k]));
      yk := exp(X[k]);
      printf("   x = %4.2f  %8.5f   %8.5f  %8.5f\n", X[k], yk, u_g(k), evalf(err[k]));
  end do:

        x      u_exact     u_g      error
   x = 0.00   1.00000    3.00000   0.00000
   x = 0.33   1.39561    7.00000   0.04883
   x = 0.67   1.94773   13.00000   0.16338
   x = 1.00   2.71828   21.00000   0.28172

  
> 

 

NULL

Download pf.mw

See...

tomleslie 13876
April 28 2022
1 0

the attached

> 

  restart;
#
# Define RSA encryption/decryption procedures
#
  RSA:= module()
              option package:
              export decrypt, encrypt:
            #
            # Decryption procedure
            #
              decrypt:= proc( n, d, mess )
                              local p:= ListTools:-Reverse
                                        ( parse~
                                          ( StringTools:-Explode
                                            ( convert
                                              ( convert
                                                ( d,
                                                  binary
                                                ),
                                                string
                                              )
                                            )
                                          )
                                        ),
                                    decryp:= Array(1..numelems(mess)),
                                    k, j, a;
                              for k from 1 to numelems(mess) do
                                  a[0]:= mess[k];
                                  for j from 1 to numelems(p) do
                                      a[j]:= a[j-1]^2 mod n;
                                  od:
                                  decryp[k]:= mul
                                              ( `if`( p[j]=1,
                                                      a[j-1],
                                                      1
                                                    ),
                                                j=1..numelems(p)
                                              ) mod n;
                              od:
                              return eval(decryp);
                        end proc;
            #
            # Encryption procedure
            #
              encrypt:= proc( n, e, mess )
                              local p:= ListTools:-Reverse
                                        ( parse~
                                          ( StringTools:-Explode
                                            ( convert
                                              ( convert
                                                ( e,
                                                  binary
                                                ),
                                                string
                                              )
                                            )
                                          )
                                        ),
                                        encryp:= Array(1..numelems(mess)),
                                        k, j, a;
                              for k from 1 to numelems(mess) do
                                  a[0]:= mess[k];
                                  for j from 1 to numelems(p) do
                                      a[j]:= a[j-1]^2 mod n;
                                  od:
                                  encryp[k]:= mul
                                              ( `if`
                                                ( p[j]=1,
                                                  a[j-1],
                                                  1
                                                ),
                                                j=1..numelems(p)
                                              ) mod n;
                              od:
                              return eval(encryp);
                        end proc;
        end module:

              
> 

  with(RSA):
#
# A test using data from
#
#  http://www.isg.rhul.ac.uk/static/msc/teaching/ic2/demo/42.htm
#
# With the supplied values for n, e, plaintext '19'
# should encrypt to '1113' and then decrypt back to
# '19'
#
# [n, e] form public key
#
  n:= 1189;
  e:= 9;
#
# Private key
#
  d:= e^(-1) mod NumberTheory:-Totient(n);
#
# Check that encryption/decryption in RSA module are at least
# plausible ie 19->1113->19
#
  encrypt( n, e, [19]);
  decrypt( n, d, %);

n := 1189

  

e := 9

  

d := 249

  

Vector(1, {(1) = 1113})

  

Array(%id = 36893488148090401484)

 (1)
> 

# Now the OP's problem
#
# [n, e] form public key
#
  n:= 187454749788911503119994043213682616233000961:
  e:= 23073697474256143563:
#
# Private key
#
  d:=e^(-1) mod NumberTheory:-Totient(n);

  mess:= Array( [ 99763088506717050716379498520357841718875740, 120165413768425112311337085191282286020515442,
                 130554852152844850464253187688342192964186591, 30831038803473309186944039044693095063005267,
                 43791870977158745785820581340119798552660187,  69969965042359538654291028155927734707714502,
                 33869965925699990903801791477929460546990002,  33844215332363231914638130818590411552396086,
                 61034870951899306892357141802983347226683758,  152680702248950827958295294572717135277913525,
                 109591116592896191360922558098912453060002989, 155489947789995346755966254731407091292570943,
                 30508537854841352888337386800131678916944679,  45657204694505066046067342403378233363404150,
                 122580283375455701674982424288406984864898071, 63190645718352016104824679398188140019415424,
                 17253386795565064409576412640690151622064852,  76432681365602448262013510061731737725777808,
                 4419855123469897669096067167819380964828479,   87913005918401414546998775379426321338719933,
                 155434516526735215629181544556450169276890711, 74695381359901114443168866905755723039288800,
                 246854149837977729662442017899408430336317,    121766227774612543867361050340041336298158533,
                 166784352080096633524033563364825999607054575, 72405993088552642775735865911360107823805504
              ]):

42699614996178234116195020050876593331093731

 (2)
> 

#
# Decrypt the encrypted message
#
  dc:= decrypt(n, d, mess);

Vector[row](26, {(1) = 123383489576056324277546637986124771145665845, (2) = 59294735682368916577199441832668436421082828, (3) = 23298112717951254177140959268401707403572298, (4) = 135306083483143972183588478788124750983888817, (5) = 140520692915900945459434750581781819734901716, (6) = 172369154102980821096292640110678268989345429, (7) = 1218178607689043481995281655552268370382698, (8) = 171139834562538131200367642187502200670759841, (9) = 103967079211641293829108632981149538396092372, (10) = 136441397210960373454903861552367586616715694, (11) = 27779966464851826101243372906874215004341021, (12) = 94772485209727728595027364529811909900964169, (13) = 137519809413815529474313088570439599935392653, (14) = 94986374324045810909928878439153557152083848, (15) = 136946575814149608803311889490599328267031261, (16) = 105093884307100664064917905369067564620244281, (17) = 90327014522086251943021419051466682748697636, (18) = 168950033358175495743577266927039472375559394, (19) = 15220313810351675164752457261149535488726071, (20) = 93747859663970601680285297658361543186789167, (21) = 83052987555357963520341436825476423137677509, (22) = 97622416778515499409147538406751417076782095, (23) = 22127691301011064624549104635398724905614759, (24) = 111762576928294948171616345725784289956370682, (25) = 186963600820326203974201460826734735260391273, (26) = 4405795800451})

 (3)
> 

#
# Re-encrypt the decrypted message, for
# comparison purposes

  ec:= encrypt(n, e, dc);

Vector[row](26, {(1) = 99763088506717050716379498520357841718875740, (2) = 120165413768425112311337085191282286020515442, (3) = 130554852152844850464253187688342192964186591, (4) = 30831038803473309186944039044693095063005267, (5) = 43791870977158745785820581340119798552660187, (6) = 69969965042359538654291028155927734707714502, (7) = 33869965925699990903801791477929460546990002, (8) = 33844215332363231914638130818590411552396086, (9) = 61034870951899306892357141802983347226683758, (10) = 152680702248950827958295294572717135277913525, (11) = 109591116592896191360922558098912453060002989, (12) = 155489947789995346755966254731407091292570943, (13) = 30508537854841352888337386800131678916944679, (14) = 45657204694505066046067342403378233363404150, (15) = 122580283375455701674982424288406984864898071, (16) = 63190645718352016104824679398188140019415424, (17) = 17253386795565064409576412640690151622064852, (18) = 76432681365602448262013510061731737725777808, (19) = 4419855123469897669096067167819380964828479, (20) = 87913005918401414546998775379426321338719933, (21) = 155434516526735215629181544556450169276890711, (22) = 74695381359901114443168866905755723039288800, (23) = 246854149837977729662442017899408430336317, (24) = 121766227774612543867361050340041336298158533, (25) = 166784352080096633524033563364825999607054575, (26) = 72405993088552642775735865911360107823805504})

 (4)
> 

#
# Check that the process of decryption followed
# by encryption returns the original (encrypted)
# message
#
  ArrayTools:-IsEqual( mess, ec);

true

 (5)
> 

 

> 

 


 

Download rsa2.mw

 

First 26 27 28 29 30 31 32 Last Page 28 of 207