acer

32373 Reputation

29 Badges

19 years, 334 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

In the example you mention, the solving was invoked by using the context-sensitive menu of operations in the left panel of the GUI.

That very fact is stated right there at the bottom of your example.

I suggest you read the corresponding Help pages (documentation). I've used a URL link above so that you can get to the on-line version of a relevent Help page.

The following example will produce a list of lists.

Each inner list has three components, for x,y, and z.

de := plot3d(y*sin(x), x = 0 .. 2*Pi, y = -2 .. 2):
temp := plottools:-transform((x,y,z)->[x,y,z])(de):
convert(plottools:-getdata(temp)[3],listlist);

As Preben has shown, your Question's second integral can be done directly in Maple 2022.2, and then further simplified. But I was unable to get that symbolic integral to compute directly in Maple 2022.1, nor any version of Maple 2021.

And your response to him mentioned only Maple 2021 in the context of that example. You then mentioned both Maple 2021 and Maple 2022 for the followup Statistics example. I would guess that means you have Maple 2022, but I'm not 100% sure.

So, for fun, here's how it can also be done in Maple 2021.2 (and onwards),

restart;

kernelopts(version);

`Maple 2021.2, X86 64 LINUX, Nov 23 2021, Build ID 1576349`

g := piecewise(z/t < 1, 0, z/t < 2, 10080-60480*z/t+156240*z^2/t^2-226800*z^3/t^3+202230*z^4/t^4-113400*z^5/t^5+39060*z^6/t^6-7560*z^7/t^7+630*z^8/t^8, 0)/t;

g := piecewise(z/t < 1, 0, z/t < 2, 10080-60480*z/t+156240*z^2/t^2-226800*z^3/t^3+202230*z^4/t^4-113400*z^5/t^5+39060*z^6/t^6-7560*z^7/t^7+630*z^8/t^8, 0)/t

H:=Int(g,t=1..3):

simplify(value(IntegrationTools:-Change(H,s=z/t,s)));

-(5/4)*piecewise(z <= 1, 0, z <= 2, -63*z^8+864*z^7-5208*z^6+18144*z^5-40446*z^4+60480*z^3-62496*z^2+48384*z-19659-8064*ln(z), z < 3, 5589-8064*ln(2), z < 6, 25248-8064*ln(2)+(7/729)*z^8-(32/81)*z^7+(1736/243)*z^6-(224/3)*z^5+(1498/3)*z^4-2240*z^3+6944*z^2-16128*z-8064*ln(3)+8064*ln(z), 6 <= z, 0)

Download mmcdara_int_2021.mw

Fwiw, your original Question's second integration example can also be forced in Maple 2015.2, by converting to Heaviside, veiling those, changing variables, computing value symbolic, and then reverting to piecewise. That kludgy approach produces some undefined discontinuities at a couple of points, which I know you don't like.  mmcdara_int_2015.mw

Those techniques involving change of variables don't seem to help in older versions for your followup Statistics example, which has an integral of a product of two piecewise, one with conditions on z/t and one on just z. This gets messed up in older versions (similar to how previous versions go awry on the two original integrals, I believe) and can erroneously produce just zero.

But perhaps you do have Maple 2022.2, which gets that followup symbolic integral as well. Or perhaps you'd be satisfied with much faster numeric integrations, as I mentioned elsewhere.

1) The forward slash dates from before modules. It's presence in Library procedure names is usually just a reflection of the directories (and subdirectories) in which the source files were stored on a Unix system. That was a convenient convention for naming internal routines within, say, old table-based packages.

Mostly the forward slash doesn't give special functionality to names. One exception is the case of the extension mechanism for kernel builtins, eg.   `print/foo` as a way of naming a procedure that would be used automatically for printing unvaluated function calls foo(...).

2) The left ticks are needed to let the parser distinguish between, say, the fraction foo/bar and the name `foo/bar`. Single left-ticks (aka name quotes) allow you to have names with nonalphanumeric characters, etc.

In contrast, the colon-dash :- is a special syntax for referring to module members. If you wrap that in name quotes then the parser sees it as just a name that happens to have some special characters in it, ie. not as a module member reference.

So these cases are actually consistent. A pair of wrapped left ticks veils the special syntactic meaning of forward slash / (ie. division), just as it would veil the special syntactic meaning of colondash :- (ie. module membership).

3) A module can be "appliable", in the sense that its name can be applied directly to arguments. For example fsolve is now a module, yet you can call  fsolve(cos(x)) .  That actually invokes 

     fsolve:-ModuleApply

on the arguments.

4) The showstat command is for rendering a display of a procedure's structure and body, including line numbers. You can also call `print` on a procedure.

You can also print (or make a string from) parts of a module. See the Help for a description of the various operands of a module. 

Note that showstat was not designed for modules and it's not clear how precisely it would display such, or why it ought to.

Note that you can call,

    showstat(fsolve::ModuleApply)

if kernelopts(opaquemodules) is `false`, and that were `true` you could use :- instead of :: .

[edit] IIRC you can also call `printf` with "%P" to specify the format, and get a plaintext, indented display of a procedure  eg.

   printf("%P",
              eval(fsolve:-ModuleApply));

The last time I checked  I could not do the same in one shot for a module. (That's inconvenient for me as module bodies can exceed the line length limit if printed after converting to unbroken string. I'd rather have nicer text pretty-printing than crude backslash linebreaks.)

If you have already executed this:

with(GroupTheory): with(DocumentTools):
lis := SearchSmallGroups('order' = 10 .. 100, 'soluble'):
S := seq(DrawSubgroupLattice(SmallGroup(g), 'highlight' = CompositionSeries(SmallGroup(g))), g in lis[10 .. 20]):

then you could do something like the following,

Tabulate(Matrix(4, 3, [S])):
That shows 3 plots per row, using the full window width by default. There are options to let you control the widths, visibility of separators, etc.

You can force the issue here by using the symbolic option of simplify and combine. Or you can use an assumption.

See the Help for some description of that symbolic option, here and here.

Sometimes the solve command will produce a "simplified" form (half your goal here).

There is a large set (nonzero measure in the RxR space of m and c) where the expression sqrt(m*c)/m and sqrt(c/m) are not equal. Hence it would hardly be an improvement if those two commands treated these two forms as equivalent, by default. I'd agree that it would be an improvement if the context-menu (context-panel) offered the choice,
   combine(..., symbolic)

restart;

g := sqrt(m*c)/m;

(m*c)^(1/2)/m

combine(simplify(g)) assuming m>0;

(c/m)^(1/2)

combine(simplify(g,symbolic),symbolic);

(c/m)^(1/2)

eq := m*omega^2 = c;

m*omega^2 = c

solve({eq, m>0}, omega);
combine(%) assuming m>0;

piecewise(m <= 0, [], 0 < m, [{omega = sqrt(c)/sqrt(m)}, {omega = -sqrt(c)/sqrt(m)}])

[{omega = (c/m)^(1/2)}, {omega = -(c/m)^(1/2)}]

Download symbolic.mw

I'm not sure that I fully understand what you mean by "improve" or "best".

You could use the mathfunc type. You could restrict that to not include the elementary functions.

In such a way you could decide up front what you want to accept as "elementary". And presumably that internal list of special functions is updated if something gets added to Maple's stock.

For example, (and this could be tweaked):

restart;

K := proc(ee)
  ormap(type,op~(0,indets(ee,function))
             minus {arccos,arccosh,arccot,arccoth,arccsc,arccsch,
                    arcsec,arcsech,arcsin,arcsinh,arctan,arctanh,
                    cos,cosh,cot,coth,csc,csch,exp,ln,sec,sech,
                    sin,sinh,tan,tanh},mathfunc);
end proc:

 

expr:=(sin(x)*EllipticF(k, y)/(BesselI(z, l)*arctan(y, x)))^WhittakerW(mu, nu, z);

(sin(x)*EllipticF(k, y)/(BesselI(z, l)*arctan(y, x)))^WhittakerW(mu, nu, z)

K(expr);

true

K( sin(x) + f(x) );

false

K( sin(x) + f(x) + BesselJ(0,x) );

true

Download specf.mw

Another approach is to store the loop index value k and the value of a in a table, and then turn that into a Matrix as a later step.

That allows you to easily handle the situation where some entries are omitted (because, say, they fail some conditional test). You can still have easy access to the accepted indices and values, and plot just those.

Here's an example. There are other ways to accomplish the same thing, some of which are terser, but I've tried to keep it understandable.

restart;

Sol:=proc(  )

  local a,k,tab;

  for k from 1 to 23 do

    a:=evalf(sin(Pi/2+9*k/23));

    if a < 0.1 then

      next;

    end if;

    tab[k]:=a;

   end do;

  return Matrix([seq([t,tab[t]],
                     t=sort([indices(tab,nolist)]))]);

end proc:

M := Sol();

M := Matrix(10, 2, {(1, 1) = 1, (1, 2) = .9244123753, (2, 1) = 2, (2, 2) = .7090764791, (3, 1) = 3, (3, 2) = .3865457698, (4, 1) = 13, (4, 2) = .3658700890, (5, 1) = 14, (5, 2) = .6931657010, (6, 1) = 15, (6, 2) = .9156718148, (7, 1) = 16, (7, 2) = .9997510142, (8, 1) = 17, (8, 2) = .9326926044, (9, 1) = 18, (9, 2) = .7246341575, (10, 1) = 19, (10, 2) = .4070289609})

(1)

plots:-display(
  plot(sin(Pi/2+9*k/23), k=1..23, color=blue),
  plot(M, style=point, symbolsize=20,
          symbol=solidcircle));

 

 

Download table_aug_Matrix.mw

 

I'm not sure how you're expecting the omega(t) values to be used. But this may be a start, or give you some ideas.

The key part here is that I'm using my earlier idea (see Reply above) of using a custom procedure F which can generate each frame.

Once you have such a procedure it is easier to test. You can call it with different arguments and see how it looks. When satisfied you can easily animate it.

[edit] I have adjusted the shape that I plotted with my method.

restart

with(plots)

r1 := .35m := 626g := -9.807
u0 := -20
v0 := -20x0 := 0z0 := 15omega0 := -24J := (2/5)*m*r1^2

eqx1 := m*(diff(u(t), t)) = 0.; eqz1 := m*(diff(v(t), t)) = m*g

eqx0 := diff(x(t), t) = u(t); eqz0 := diff(z(t), t) = v(t)

eqr1 := J*(diff(omega(t), t)) = 5000*sin(t)

ini := x(0) = x0, u(0) = u0, z(0) = z0, v(0) = v0, omega(0) = omega0

sol1 := dsolve({eqr1, eqx0, eqx1, eqz0, eqz1, ini}, {omega(t), u(t), v(t), x(t), z(t)}, type = numeric, output = listprocedure)

animx := subs(sol1, x(t)); animz := subs(sol1, z(t)); animomega := subs(sol1, omega(t))

F := proc (x, z, r, q) plots:-display(plottools:-rotate(plots:-display(plottools:-sector([x, z], q, 0 .. Pi, color = "LightBlue"), plottools:-sector([x, z], q, Pi .. 2*Pi, color = grey)), r, [x, z])) end proc

animate(F, [animx(t), animz(t), animomega(t), .5], t = 0 .. 1, scaling = constrained, frames = 150)

Download how_rot_ac.mw

By calling parse you turn the strings back into Maple expressions, so naturally when you print them as normal output then you get the same effect as for the original, ie. sequences are pretty-printed with spaces after the commas, whether they are in a list or just as a sequence.

If you want to line-print strings, without the quotation marks appearing, then you could use printf (on the strings, with %s format) instead of print or lprint. Note that conversion to strings might be unnecessary -- your particular goal is still unclear.

If you want to export the entries to an external file then you could use the ExportMatrix command, eg.,
    ExportMatrix("myfile.csv", Matrix([ans1]));
That produces a file named myfile.csv with three lines, each with the entries separated by commas and without spaces. That command has some additional flexibility along those aspects.

What are you actually trying to accomplish with sequences of numbers displayed without spaces after the commas? Are you trying to export it to an external file, for some other software? Or are you going to copy&paste from what is displayed on your Maple GUI?

It is a very common phenomenon on this forum for people to have some task they want to solve, then ask about only part of their approach to solving it, without describing the actual underlying goal or full process that they want. This usually leads to guesswork and unnecessary and inefficient back-and-forth.

In the Maple programming language (1D plaintext, also known as Maple Notation) the word for infinity is:

   infinity

That entry box in the popup is a Maplets TextField, and it expects plaintext Maple Notation.

This seems to be a consequence of the result from DiagonalMatrix having shape=diagonal and storage=diagonal.

A workaround is to force off the compact shape/storage (see two ways, below).

I will submit a bug report.

restart

with(LinearAlgebra); alias(`&bigotimes;` = LinearAlgebra:-KroneckerProduct); interface(rtablesize = 16); kernelopts(version)

`&bigotimes;`

[10, 10]

`Maple 2021.2, X86 64 LINUX, Nov 23 2021, Build ID 1576349`

Matrices nonconformable, Maple should give error message:

`&bigotimes;`(Matrix(1, 4, 1), IdentityMatrix(4)); Dimension(%); M := Matrix(4, 4, shape = symmetric, symbol = m); Dimension(%)

Matrix([[1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0], [0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0], [0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0], [0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1]])

4, 16

M := Matrix(4, 4, {(1, 1) = m[1, 1], (1, 2) = m[1, 2], (1, 3) = m[1, 3], (1, 4) = m[1, 4], (2, 2) = m[2, 2], (2, 3) = m[2, 3], (2, 4) = m[2, 4], (3, 3) = m[3, 3], (3, 4) = m[3, 4], (4, 4) = m[4, 4]}, storage = triangular[upper], shape = [symmetric])

4, 4

MatrixMatrixMultiply(`&bigotimes;`(Matrix(1, 4, 1), IdentityMatrix(4)), DiagonalMatrix(Diagonal(M), shape = []))

Error, (in LinearAlgebra:-Multiply) first matrix column dimension (16) <> second matrix row dimension (4)

NULLMatrixMatrixMultiply(`&bigotimes;`(Matrix(1, 4, 1), IdentityMatrix(4)), DiagonalMatrix(Diagonal(M), storage = rectangular))

Error, (in LinearAlgebra:-Multiply) first matrix column dimension (16) <> second matrix row dimension (4)

`&bigotimes;`(Matrix(1, 4, 1), IdentityMatrix(4)).DiagonalMatrix(Diagonal(M))

Matrix([[m[1, 1], 0, 0, 0], [0, m[2, 2], 0, 0], [0, 0, m[3, 3], 0], [0, 0, 0, m[4, 4]]])

Download Starnge_ac.mw

Is you lprint the SearchSet you can see that its float values are not what appears in the 2D Input.

I don't know how you managed to get that effect (cut&paste, something with numeric formatting?).

You could re-enter it manually.

But you may still be faced by discrepencies if you expect floats to be matched in the case that they have differing numbers of trailing zeroes.

It might be better if you described precisely what kind of matching you are after, ie. float comparisons to a fixed precision, etc. It may be that verify it a more suitable tool for the comparison, but that depends on what comparison you are after.

Also, you asked about why some sorting occurs. That's how Maple acts with sets of floats. If you don't want them sorted then you could use lists of Vectors instead.

restart; with(ListTools)

Pent := [Vector(3, {(1) = 2.0646, (2) = -1.5000, (3) = 3.3405}), Vector(3, {(1) = -.78860, (2) = -2.4271, (3) = 3.3405}), Vector(3, {(1) = -2.5520, (2) = 0.8100500000e-10, (3) = 3.3405}), Vector(3, {(1) = -.78860, (2) = 2.4271, (3) = 3.3405}), Vector(3, {(1) = 2.0646, (2) = 1.5000, (3) = 3.3405})]

[Vector[column](%id = 36893628814707266132), Vector[column](%id = 36893628814707266252), Vector[column](%id = 36893628814707266372), Vector[column](%id = 36893628814707266492), Vector[column](%id = 36893628814707266612)]

PentSet := ([seq])(convert(Pent[i], set), i = 1 .. 5)

[{-1.5000, 2.0646, 3.3405}, {-2.4271, -.78860, 3.3405}, {-2.5520, 0.8100500000e-10, 3.3405}, {-.78860, 2.4271, 3.3405}, {1.5000, 2.0646, 3.3405}]

PentSetSort := ([seq])(sort(PentSet[i]), i = 1 .. 5)

[{-1.5000, 2.0646, 3.3405}, {-2.4271, -.78860, 3.3405}, {-2.5520, 0.8100500000e-10, 3.3405}, {-.78860, 2.4271, 3.3405}, {1.5000, 2.0646, 3.3405}]

SearchSet := {-2.427050982, -.7885966699, 3.340549093}; BinarySearch(PentSetSort, SearchSet, `subset`); BinarySearch([{-1.500000000, 2.064572880, 3.340549093}, {-2.427050982, -.7885966699, 3.340549093}, {-2.551952424, 3.340549093, 8.100498050*10^(-11)}, {-.7885966673, 2.427050983, 3.340549093}, {1.500000001, 2.064572880, 3.340549093}], {-2.427050982, -.7885966699, 3.340549093}, `subset`)

0

2

lprint(SearchSet)

{-2.427050982, -.7885966699, 3.340549093}

SearchSet := {-2.4271, -.7886, 3.3405}

{-2.4271, -.7886, 3.3405}

BinarySearch(PentSetSort, SearchSet, `subset`)

0

SearchSet := {-2.4271, -.78860, 3.3405}

{-2.4271, -.78860, 3.3405}

BinarySearch(PentSetSort, SearchSet, `subset`)

2

lprint(PentSetSort)

[{-1.5000, 2.0646, 3.3405}, {-2.4271, -.78860, 3.3405}, {-2.5520, .8100500000e-\
10, 3.3405}, {-.78860, 2.4271, 3.3405}, {1.5000, 2.0646, 3.3405}]

([seq])(convert(Pent[i], set), i = 1 .. 5)

[{-1.5000, 2.0646, 3.3405}, {-2.4271, -.78860, 3.3405}, {-2.5520, 0.8100500000e-10, 3.3405}, {-.78860, 2.4271, 3.3405}, {1.5000, 2.0646, 3.3405}]

([seq])(convert(Pent[i], list), i = 1 .. 5)

[[2.0646, -1.5000, 3.3405], [-.78860, -2.4271, 3.3405], [-2.5520, 0.8100500000e-10, 3.3405], [-.78860, 2.4271, 3.3405], [2.0646, 1.5000, 3.3405]]

 

Download Binary_search_ac.mw

An explicitl description of what you're actually trying to accomplish would be a good start to getting best advice.

It would be more helpful and direct if you were to supply the complete example when first submitting the Question.

Hopefully this uses the dsolve solution to produce your line animation in a satisfactory manner.

restart

q1 := 5b1 := 18l1 := 16r1 := .35ms := 626g1 := -9.807
u0 := -19.4

v0 := -17.7x0 := 0z0 := 15xa := 0za := z0+1mu1 := .35omega0 := -24

mn := b1*g1*q1*za+g1*mpfz := mn/zaJs := (2/5)*ms*r1^2

"Lp(t):=sqrt((xn(t)-xa)^(2)+(zn(t)-za)^(2)) :"

"Lr(t):=piecewise((za-Lp(t))>0,za-Lp(t),0) :"

"A1(t):=piecewise((xa-xn(t))>0,xa-xn(t),0) :"

"sinalpha1(t):=piecewise((xa-xn(t))>0,(xa-xn(t))/(Lp(t)),0) :"

"cosalpha1(t):=piecewise((xa-xn(t))>0,((za-zn(t)))/(Lp(t)),1) :"

"alpha2(t):=arccos(cosalpha1(t))+phi(t) :"

``

"S1(t):=piecewise(zn(t)>=0 and xn(t)<=xa,fz*Lr(t),0):"

"S1z(t):=S1(t)*cos(phi(t)):"

"S1x(t):=S1(t)*sin(-phi(t)):"

"S2(t):=piecewise(xn(t)<=xa and zn(t)>0, (-1) S1(t)*exp(mu1*alpha2(t)),0):"

"S2x(t):=S2(t)*sinalpha1(t):"

"S2z(t):=S2(t)*cosalpha1(t):"

"Sx(t):=S2x(t)+S1x(t):"

"Sz(t):=S1z(t)+S2z(t)+ms*g1:"

shape_factor := 0

"Fr(t):=S2(t)+S1(t)+shape_factor:"

NULL

"m(t):=Lr(t)q1*b1:"

mp := 300

"Jp(t):=1/(3)*m(t)*Lr(t)^(2)+mp*Lr(t)^(2) :"

"M(t):=g1*m(t)*(Lr(t))/(2)*sin(phi(t))+g1*mp*sin(phi(t)):"

k := 800*ms

NULL

"Es(t):=(ms*g1*(z0-zn(t)))+1/(2)*ms*(vn(t)^(2)+un(t)^(2))+1/(2)*Js*omega(t)^(2):"

"En(t):=(1)/(2)*Jp(t)*((&DifferentialD;)/(&DifferentialD; t)phi(t))^(2)+-(m(t)+mp)*g1*(Lr(t)-cos(phi(t))*Lr(t)):"

eqx1 := ms*(diff(u(t), t)) = 0.; eqz1 := ms*(diff(v(t), t)) = ms*g1

eqx0 := diff(x(t), t) = u(t); eqz0 := diff(z(t), t) = v(t)

``

eqx1n := ms*(diff(un(t), t)) = Sx(t); eqz1n := ms*(diff(vn(t), t)) = ms*g1+S1z(t)+S2z(t)

eqx0n := diff(xn(t), t) = un(t); eqz0n := diff(zn(t), t) = vn(t)

eqr1 := Js*(diff(omega(t), t)) = (S2(t)+S1(t)+shape_factor)*r1

``

u0n := u0*ms/(ms+(1/3)*m(t)+mp)

NULL

u0n1 := u0*ms/(ms+mn)

phid00 := u0n1/(l1-z0)

phid0 := .7*eval[recurse](-u0n(t)/(l1-z0), [t = 0, xn(0) = x0, zn(0) = z0])

simplify(phid0)

"znb(t):=(Lr(t)-(sin((Pi)/(2)- phi(t))*Lr(t))):"

"xnb(t):=Lr(t)*cos(Pi/(2)-phi(t)):"

eqp := (diff(phi(t), t, t))*Jp(t) = M(t)-k*(diff(phi(t), t))

with(plots)

ini := x(0) = x0, u(0) = u0, z(0) = z0, v(0) = v0, xn(0) = x0, un(0) = u0, zn(0) = z0, vn(0) = v0, omega(0) = omega0, phi(0) = 0, (D(phi))(0) = phid0

sol1 := dsolve({eqp, eqr1, eqx0, eqx0n, eqx1, eqx1n, eqz0, eqz0n, eqz1, eqz1n, ini}, {omega(t), phi(t), u(t), un(t), v(t), vn(t), x(t), xn(t), z(t), zn(t)}, type = numeric, output = listprocedure)

XNB := unapply(subs(xn = XN, zn = ZN, phi = PHI, xnb(t)), t)

ZNB := unapply(subs(xn = XN, zn = ZN, phi = PHI, znb(t)), t)

XN := eval(xn(t), sol1)

ZN := eval(zn(t), sol1)

PHI := eval(phi(t), sol1)

NULL

NULL

NULL

NULL

plots:-animate(plots:-display, [plottools:-line([xa, za], [XN(t), ZN(t)], color = red, thickness = 3), plottools:-line([XN(t), ZN(t)], [XNB(t), ZNB(t)], color = blue, thickness = 3), plot([[XN(t), ZN(t)]], style = point, symbol = solidcircle, symbolsize = 15, color = red), plot([[XNB(t), ZNB(t)]], style = point, symbol = solidcircle, symbolsize = 15, color = blue)], t = 0 .. .7, background = plots:-pointplot([[xa, za]], symbol = solidcircle, symbolsize = 15, color = green))

 

Download anim_ac.mw

You have b assigned as a Matrix. That's your basic mistake, since you were only indexing into it as b[i], as if it had just one dimension.

Try making it b as a Vector instead, ie,

   b := Vector(3, [-1, 7, -7]);

Below, I have removed also the loading of the VectorCalculus package since that rebinds the Vector command which is hardly useful here. I changed the call of Norm to be explicitly a call to LinearAlgebra:-Norm.

I suggest that you check that it works as you intended.

restart

with(plots):

``

SOR := proc (n, A, b, x, w, max_iter, tol) local i, j, k, xGS, err, x1; err := 1; while tol < err do x1 := x; for k to max_iter do for i to n do xGS := (b[i]-add(Typesetting:-delayDotProduct(A[i, j], x[j]), j = 1 .. i-1)-add(Typesetting:-delayDotProduct(A[i, j], x[j]), j = i+1 .. n))/A[i, i]; x[i] := (1-w)*x[i]+w*xGS end do end do; err := LinearAlgebra:-Norm(x-x1); x1 := x end do end proc:

``

A := Matrix(3, 3, [[3, -1, 1], [-1, 3, -1], [1, -1, 3]])

A := Matrix(3, 3, {(1, 1) = 3, (1, 2) = -1, (1, 3) = 1, (2, 1) = -1, (2, 2) = 3, (2, 3) = -1, (3, 1) = 1, (3, 2) = -1, (3, 3) = 3})

(1)

b := Vector(3, [-1, 7, -7]); x := `<,>`(1, 1, 1); x := convert(x, Vector)

b := Vector(3, {(1) = -1, (2) = 7, (3) = -7})

 

x := Vector(3, {(1) = 1, (2) = 1, (3) = 1})

(2)

n := 3;

3

 

.9

 

7

 

0.1e-3

(3)

ans := SOR(3, A, b, x, .9, 15, 0.1e-3)

ans := Vector(3, {(1) = .9999999896, (2) = 1.999999980, (3) = -2.000000006})

(4)

A.ans-b

Vector(3, {(1) = -0.1719999987e-7, (2) = -0.4360000005e-7, (3) = -0.8399999807e-8})

(5)

NULL

Download code_zeineb_ac.mw

First 59 60 61 62 63 64 65 Last Page 61 of 336