@Carl Love

That's crazy. Thanks for the info. Where did you read this? As of right now I have changed some things so that the code works for all types of equations. It looks like this if you are interested.

restart:

Digits:=12:

Runtime:=time():

with(LinearAlgebra):

with(plots):

with(orthopoly):

with(ArrayTools):

N := 2:

phix := Vector(N):

X := Vector(N):

phix[1] := cos(x[2]):

phix[2] := 3*cos(x[1]):

X[1] := 1.0:

X[2] := 2.0:

#Code AndersonAcceleration.

AndersonAcceleration:=proc(N,phi,X)

global x, xS, here1, here2;

local mMax, itmax, atol, rtol, droptol, beta, AAstart, res_hist, df, DGg, gamma, X0;

local DG, mAA, iter, gval,fval,res_norm, condDF, tol, f_old, g_old, y, i, k, Q, R, QRg, dfT, DGT;

(*

%------------------------------------------------------------------------

% Function xS = AndersonAcceleration(N,phi,x0).

%

% Fixed-point iteration with or without Anderson acceleration.

% 'phi' is the fixed point iteration map and 'xS' is the

% solution, so that xS = phi(xS).

%

% Input arguments:

% X0 = Initial value solution. Note: In this function this variable

% must be a column vector.

% 1. 'mMax' = maximum number of stored residuals (non-negative integer).

% NOTE: 'mMax' = 0 => no acceleration. default=1000

% 2. 'itmax' = maximum allowable number of iterations. default=1000

% 3. 'atol' = absolute error tolerance. default=1.0e-6

% 4. 'rtol' = relative error tolerance. default=1.0e-3

% 5. 'droptol' = tolerance for dropping stored residual vectors to

% improve conditioning: If 'droptol' > 0, drop residuals if the

% condition number exceeds droptol; if droptol <= 0,

% do not drop residuals.

% 6. 'beta' = damping factor: If 'beta' > 0 (and beta ~= 1), then

% the step is damped by beta; otherwise, the step is not damped.

% NOTE: 'beta' can be a function handle; form beta(iter), where iter

% is the iteration number and 0 < beta(iter) <= 1.

% 7. 'AAstart' = acceleration delay factor: If 'AAstart' > 0, start

% acceleration when iter = AAstart.

%

% Output:

% xS = Solution vector.

%

% The notation used is that of H.F. Walker: Anderson Acceleration:

% Algorithms and implementation

%------------------------------------------------------------------------

*)

mMax := 100:

itmax := 100:

atol := 1.0e-8:

rtol := 1.0e-6:

droptol := 1.0e12:

beta := 1.0:

AAstart := 0:

# Initialize storage arrays and number of stored residuals.

DG := Matrix():

df := Matrix():

DGg := Vector(N);

QRg := Vector(N);

mAA := 0:

X0 := X;

for iter from 0 to itmax do

x:=X0:

gval := Vector(phi):

fval := gval - X0:

res_norm := norm(fval,2):

# Set the residual tolerance on the initial iteration.

if iter = 0 then

tol := max(atol,rtol*res_norm):

fi:

# Convergence test, if converged the loop stops.

if res_norm <= tol then

print('iter',iter,res_norm);

break; # Breaks for-loop

fi:

# If resnorm is larger than 1e8 at iter > 5, problem stops

if res_norm >1e8 and iter > 5 then

print('no_convergence',res_norm);

break; # Breaks for-loop, diverged

fi:

# Fixed point iteration without acceleration, if mMax == 0.

if mMax = 0 or iter < AAstart then

# We update E <- g(E) to obtain the next approximate solution.

X0 := gval:

else

# With Anderson acceleration.

# Update the df vector and the DG array.

if iter > AAstart then

if mAA < mMax or Size(df,2) = 1 then

df := Concatenate(2,df,fval-f_old):

DG := Concatenate(2,DG,gval-g_old):

else

df := Concatenate(2,df[..,-1],fval-f_old):

DG := Concatenate(2,DG[..,-1],gval-g_old):

fi:

mAA := mAA + 1:

fi: # iter

# We define the old g and f values for the next iteration

f_old := fval;

g_old := gval;

if mAA = 0 then

# Initialization

# If mAA == 0, update X <- g(X) to obtain themah next approximate

# solution. No least-squares problem is solved for iter = 0

X0 := gval;

else

if mAA > 1 then

Q,R := QRDecomposition(df,datatype=float,output='NAG');

if type(Q, 'Matrix'(square)) then

condDF := ConditionNumber(Q);

else

condDF := 1;

fi:

while condDF > droptol and mAA > 1 do

if mAA = 2 then

df := convert(df[..,2..-1],Vector);

DG := convert(DG[..,2..-1],Vector);

else

df := df[..,2..-1];

DG := DG[..,2..-1];

fi:

Q,R := QRDecomposition(df,datatype=float,output='NAG');

mAA := mAA - 1;

if type(Q, 'Matrix'(square)) then

condDF := ConditionNumber(Q);

else

condDF := 1;

fi:

od:

if Size(df,2) > 1 then

gamma := LeastSquares([Q,R],fval);

else

R := norm(df,2);

Q := R/~df;

gamma := R/~Transpose(Q).fval;

fi:

else

R := norm(df,2);

Q := R/~df;

gamma := R/~Transpose(Q).fval;

fi:

if Size(gamma,1) > 1 then

DGg:=DG.gamma:

else

DGg:=DG*gamma;

fi:

# Update the approximate solution.

X0 := gval - DGg;

# Damping for non-zero beta

if beta > 0 and beta <> 1 then

if mAA = 1 then

QRg := Q*R*gamma;

else

QRg := df.gamma;

fi:

X0 := X0 - (1-beta)*(fval - QRg);

fi:# isa(beta ...

fi: # mAA = 0

fi:# mMax == 0

od:

X[1..N] := X0[1..N];

end:

AndersonAcceleration(N,phix,X):