## 8333 Reputation

11 years, 118 days

## The attached code works...

in Maple18 animMaple18.mw

in Maple 2015 animMaple2015.mw

in Maple 2016 animMaple2016.mw

in Maple 2017 animMaple2017.mw

in Maple 2018 animMaple2018.mw

in Maple 2019 animMaple2019.mw

in Maple 2020 animMaple2020.mw

in MAple 2021 animMaple2021.mw

I'm afraid that I have no conveneinent way of checking the function of this code  in Mpale versions older than Maple 18 (released in 2014). I only have the last eight versions of Maple loaded on this machine. Any version more than eight years old, I can't really help with (other than advise you to update!)

## Hmmm...

This site is about using Maple to solve problems so first question: what exactly is Maple doing incorrectly?

If I go through your previous response more-or-less line by line, then you have to ask yourself some pretty basic questions, which I have itemised in the following

The general solver -LinearSolve() confirms the existence of a solution.

Incorrect. statement. LinearAlgebra:-LinearSolve demonstrate that there are an infinite number of solutions

The -IterativeFormula() suggests that this solution can be reached using one of the standard iterative methods (such Gauss-Seidel, Jacobi, etc).

Incorrect statement. Since there are an inifinite number of solutions, all the IterativeFormula() command does is find one of them

I guess, the fact that these methods can deal with the case of spectral radius of M being 1 is not so surprizing, since the singular matrix (I-M) from equation (5) of my write-up can be inverted using the Moore-Penrose pseudoinverse and still generate an (up-to-scale) unique solution.

Overcomplicating the issue: For a numerical process to generate one (of an infinite number of) solution(s), all it has to do is to assign an arbitrary value to one of the variables, throw away one of the redundant equations, and solve the remaining three independent (consistent) linear equations in three unknowns. Not exactly difficult

Presumably, however, the -IterativeFormula() applies a different type of iteration then the one you've coded in Example1.mw and Example2.mw above (where we start with an initial, model-consistent guess for the unknowns (of 0's) on the right hand side (RHS), calculate the left-hand side (LHS), plug those values again on the RHS, etc.), for which we have a convergence in Ex1 and a divergence in Ex2. (Presumably, the -IterativeFormula() plugs in the LHS values with some sort of 'damping' factor on the RHS to ensure the convergence to the existent solution. However, the existence of a solution that can be reached iteratively, e.g., using -IterativeFormula(), does not yet imply the stability of the equilibrium, right?).

Complete misrepresentation: the problem you posted in your very first worksheet was for a system of non-linear equations - lots of products of y[1]*y[2] etc. I was doubtful about the success of an iterative process in solving this non-linear system, which is why my original answer was headlined "use with care". At some point in the evolution of this thread, you morphed from a system of non-linear equations to a system of linear equations - I have no idea how/why this occurred, but if your original post had specified a system of linear equations, then no way would I have proposed the iterative solution in my original response. For systems of linear equations, there are several "better" methods of constructing an iterative solution - see for example the IterativeFormula() command

Hence, I wonder whether the fact of a divergence in Example2.mw (via a "mechanical" updating of the RHS using LHS) is some sort of indication that the equilibrium in Example2 is unstable and, if so, how can I formally tackle this question?

More-or-less meaningless: As stated above, I have little confidence that the iterative approach in my original response would work for non-linear system. For linear systems, far better iteration processes are avaliable. (See the documentation for the IterativeFormula() command.) I also have no idea what you mean by terms such as "equilbrium" or "unstable" in this context. All you have is two linear systems of four non-independent equations, each of which has an infinite number of solutions

## Interesting...

As I pointed out in an earlier response, bith ot the examples you quote can be solved using the

LInearAlgebra:-LinearSolve()

command, providding (in both cases) an infinite family of solutions because (for both examples) the coefficient matrix is of rank=3 and the systems have four unknowns.

One can also examine iterative solutions of both examples using the

Student[NumericalAnalysis]:-IterativeFormula()

command. One interesting output of this command is that it returns the spectral radius for both systems as identically 1 which suggests(?) that an iterative process ought not to work(?). However in both cases if one examines the iterates, then a solution is obtained for both examples - and in very few iterations.

I can only suggest that you "play around" a little with IterativeFormula() command (there a lot of options to try) in order to see if this provides you with any enlightenment.

See the attached

  restart;
ex1 := {y[1] = 109/2000 + (327*y[4])/400 + (507*y[2])/400 - (327*y[3])/400 - (107*y[1])/400,
y[2] = -109/3000 - (109*y[4])/200 + (31*y[2])/200 + (109*y[3])/200 + (169*y[1])/200,
y[3] = (3*y[1])/10 + y[2]/5 - 1/50 + (3*y[3])/10 + y[4]/5,
y[4] = y[1]/10 + (2*y[2])/5 - 1/150 + y[3]/10 + (2*y[4])/5}:
ex2 := {y[1] = (109*y[4])/100 + (77*y[2])/50 - (109*y[3])/100 - (27*y[1])/50 + 109/2000,
y[2] = -(109*y[4])/150 - (2*y[2])/75 + (109*y[3])/150 + (77*y[1])/75 - 109/3000,
y[3] = (3*y[1])/10 + y[2]/5 - 3/200 + (3*y[3])/10 + y[4]/5,
y[4] = y[1]/10 + (2*y[2])/5 - 1/200 + y[3]/10 + (2*y[4])/5}:
#
# Solve the linear system ex1 in the most "general" way
# possible, using LinearAlgebra:-LinearSolve. This
# will produce an infinite family of solutions
#
LASOL1:= evalf
( LinearAlgebra:-LinearSolve
( LinearAlgebra:-GenerateMatrix
( ex1,
union(indets~(ex1)[])
)
)
);
#
# Examine numeric solutions using the Student[NumericalAnalysis]
# package. Interestingly it report the spectral radius as 1,
# which ought(?) to mean that the iterative process should fail
#
Student[NumericalAnalysis]:-IterativeFormula
( LinearAlgebra:-GenerateMatrix
( ex1,
union(indets~(ex1)[])
),
);
#
# But examine the iterates anyway - they appear to converge
#
#
Student[NumericalAnalysis]:-IterativeFormula
( LinearAlgebra:-GenerateMatrix
( ex1,
union(indets~(ex1)[])
),
iterations=10,
initialapprox=Vector([0.,0.,0.,0.]),
output=iterates
);
#
# Verify that the (single) iterative solution returned by the
# Student[NumericalAnalysis]:-IterativeFormula command belongs
# to family of solutions returned by the LinearAlgebra:-LinearSolve()
# command
#
eval( LASOL1,
isolate
(  LASOL1[2]=0,
indets(LASOL1[2])[])
);


## More observations...

If I convert your examples to matrix form and use LinearSolve(), both produce an infinite family of solutions.

Double-checking: in both cases the rank of the coefficient matrix is 3 and the rank of the augmented matrix is also 3.

By the theorem here

https://en.wikipedia.org/wiki/Rouch%C3%A9%E2%80%93Capelli_theorem,

since the augemnted matrix and the oefficient matrix have the same rank, then the system is consistent. However since the 'rank' is less than the number of variables, no unique solution exists in either case,

 restart;
with(LinearAlgebra):
ex1 := {y[1] = 109/2000 + (327*y[4])/400 + (507*y[2])/400 - (327*y[3])/400 - (107*y[1])/400,
y[2] = -109/3000 - (109*y[4])/200 + (31*y[2])/200 + (109*y[3])/200 + (169*y[1])/200,
y[3] = (3*y[1])/10 + y[2]/5 - 1/50 + (3*y[3])/10 + y[4]/5,
y[4] = y[1]/10 + (2*y[2])/5 - 1/150 + y[3]/10 + (2*y[4])/5}:
ex2 := {y[1] = (109*y[4])/100 + (77*y[2])/50 - (109*y[3])/100 - (27*y[1])/50 + 109/2000,
y[2] = -(109*y[4])/150 - (2*y[2])/75 + (109*y[3])/150 + (77*y[1])/75 - 109/3000,
y[3] = (3*y[1])/10 + y[2]/5 - 3/200 + (3*y[3])/10 + y[4]/5,
y[4] = y[1]/10 + (2*y[2])/5 - 1/200 + y[3]/10 + (2*y[4])/5}:
LinearSolve(GenerateMatrix( ex1, union(indets~(ex1)[]) ));
LinearSolve(GenerateMatrix( ex2, union(indets~(ex2)[]) ));

(A,b):=GenerateMatrix(ex1,union(indets~(ex1)[]) ):
Rank(A);
M:=GenerateMatrix( ex1, union(indets~(ex1)[]), augmented=true ):
Rank(M);
(A,b):=GenerateMatrix(ex2, union(indets~(ex2)[])):
Rank(A);
M:=GenerateMatrix( ex1, union(indets~(ex2)[]), augmented=true ):
Rank(M);


See the attached rank.mw

## Some observations...

If I use the inbuilt 'fsolve()' command and specify the "starting value", then answers are obtained for both of your examples. For the first example, the answer is

{ y[1] = 0.05138004246 + 1.000000000*y[4],
y[2] = 1.000000000*y[4] + 0.005095541402,
y[3] = -0.005095541405 + 1.000000000*y[4],
y[4] = y[4]}

So y[4] is arbitrary, and other variables can be determined in terms of y[4]. The same answer is obtained for a couple of different starting values. Whether this answer would be obtained for all starting values is a much more difficult question. See the worksheet ex1.mw

For the second example 'fsolve' obtains the answer

{y[1] = -2.174290792, y[2] = -2.211876999, y[3] = -2.218083895, y[4] = -2.214980447}

Again the same answer is obtained from a couple of different starting values. See the worksheet ex2.mw

## OOooops...

@vv Of course you are correct!

I can't really work out what this worksheet is trying to achieve - but I think you have (inadvertently) stumbled across a bug in Maple. This potential bug is so serious that I have posted it as a separate question on this site

seeking any comments, before I submit an SCR to MapleSoft

For your current problem I strongly advise that you do not use indexed variable names (ie X[2], X[4]) as the 'loop index' in any seq(), mul() or add() command. There appears to be no reason for using indexed variables in your worksheeet, so just change them to 'atomic' variables (ie X__2, X__4 etc)

## You state...

I want to solve the following system (please see the code)

What code??

Upload a worksheet using the big green up-arrow in the Mapleprimes toolbar

## Idle curiosity...

used the DOS command line, and turned off the warning messages using DOS option, and now it finished almost instantly.

Obvious question - what happens if you turn off warning messages in the GUI?

## Maybe...

the attached is what you want. Note that I have set the plot ranges and view ranges, purely for illustrative purposes - change these to anything you want.

 > restart;   gamma__1:= 2:   expr:=gamma__1*(100-x)/100-x*0.05*(gamma__1*(100-x)/100):   k:= unapply       ( piecewise         ( or( entries~( solve(expr>0, [x]),                           'nolist'                         )[]               ),           expr,           0         ),         x       );   M:= 41539.42878*( 0.010+0.0001*((k(x)+2*gamma__1)/2)^2 ):   plot( k,         -100..200,         view=[-100..200, -5..25],         axes=boxed       );   plot( abs(M),         x=-100..200,         view=[-100..200, 0..1200],         axes=boxed       );
 >

## Of course you could use...

a small variation on my original suggestion. I regret I did not handle the possibility of having to deal with "escape" characters, ie "\"  in my original response

However for the text file

1,2,"0",4
1,2,"1",4
1,2,"A",4
1,2,"\int \frac{(a+b x)^2 (A+B x)}{\sqrt{c+d x} \sqrt{e+f x} \sqrt{g+h x}} \, dx",4

the attached worksheet would seem to provide the "correct" answer

 > restart:   fName:="C:/Users/TomLeslie/Desktop/data2.txt":   myRead:= proc( myFile )                  local j:=1,                        l,                        line;                  uses StringTools;                  while true do                        line:=readline( myFile );                        if   line <> 0                        then l[j]:=[parse(SubstituteAll(line, "\\", "\\\\"))];                             j:=j+1;                        else break;                        fi;                 od:                 return Matrix([entries(l, 'nolist')])          end proc:    M:=myRead( fName )
 (1)
 >

In future please try to avoid asking a question - then pointing out that perfectly adeuate answers are not "appropriate" simply because you were incapable of formulating the question correctly

## Makes no sense...

You have defined the function L() as

L:=z->piecewise( d<=z,
1-2*xi*(cos((2*3.14)*((z-d)*(1/2))-1/4)-(7/100)*cos((32*3.14)*(z-d-1/2))),
z<=d+1,
1-2*xi*(cos((2*3.14)*((z-d)*(1/2))-1/4)-(7/100)*cos((32*3.14)*(z-d-1/2))),
1
);

So one can compute L(0.71) - the answer is 0.8074456472.

When you say that you want

L(0.71)=0

then you are saying 0.8074456472=0 which makes no sense

## Well...

This latets piece of code is non-executable - maybe you think people here have nothing better to do than stitch together parts of yur first worksheet with parts of your second worksheet, just in order to get something which will even execute. The bad news is that this *probably* isn't going to happen. It is your responsibility to provide a complete executbale worksheet whihc illustrates your problem, prefereably using the big green up-arrow in the Mapleprimes toolbar.

If I make some random guesses about the information required by your second worksheet, then the solve() command appears to always generate two solutions.

Why do you think that these will be in some particular order?

What order do you want them in?

How is Maple supposed to "know" what order you want them in? Maybe you think that MAple can guess what you want?

## More like...

the attached?

 > eqns:=[2*x+3*cos(y)+z^2=0, sin(x)+y+2*z=1];
 (1)
 > plot( [ 'eval'(x,'fsolve'(eval(eqns, z=Z))),         'eval'(y,'fsolve'(eval(eqns, z=Z)))       ],       Z=-5...5,       color=[red,blue],       thickness=4     );
 >