Iterative method for fixed point of a mapping with domain and range contained in $\mathbb{R}^{2}$

Let $$X$$ be a metric space and $$T:X\rightarrow X$$ be a function. A point $$x\in X$$ is called a fixed point of $$T$$ if $$x=Tx$$. For example, if $$X=[0,1]$$ and $$Tx=\frac{x}{2}$$ then $$T$$ has a unique fixed point $$x=0$$. Now I can compute this fixed by the iterative method $$x_{n}=Tx_{n-1}$$ where $$n=1,2,3,...$$ and $$x_{0}\in X$$ is any initial guess. My code is the following

T[x_] := T[x] = (x/2)
x[0] = 0.1;
x[n_] := x[n] = T[x[n - 1]]
NumberForm[a1 = Table[x[i], {i, 0, 20}], 4]


After run this, I get the values as

0.1, 0.05, 0.025, 0.0125, 0.00625, 0.003125, 0.0015625


Now I have a problem in coding and no idea how I set the code similar to my above code when the domain and range are in $$\mathbb{R}^{2}$$. For example if $$X=[0,1]\times [0,1]$$ and set $$T:X\rightarrow X$$ by $$T(x,y)=(\frac{x}{2},\frac{y}{3})$$ then $$(0,0)$$ is the unique fixed point of $$T$$. Now how is compute this fixed point by the iterative method above in this case in mathematica?

• You can represent points as lists with 2 elements, e.g. T[{x_,y_}] := {x/2,y/3}. Commented Nov 19, 2022 at 10:18

• Iterative the mapping Apply[{x, y} |-> {x/2, y/3}].
pts = NestList[Apply[{x, y} |-> {x/2, y/3}], {20., 10.}, 10]
Arrow /@ Partition[pts, 2, 1]}, Axes -> True]


• Calculate the fixed point by DiscrteLimit.
DiscreteLimit[
RSolveValue[{x[n + 1] == x[n]/2, y[n + 1] == y[n]/3, x[0] == x0,
y[0] == y0}, {x[n], y[n]}, n], n -> ∞]


{0,0}

Try

FixedPoint[(#/2) &, x, 10]
(* x/1024 *)


and play with its last argument 10, 15, 100...

Have fun!

2-D Fixed-point iteration to solve Lame-Embden integral equation with matrix multiplication

You can modify the attached code for your or another problem. The code shows how to deal with finding a root via the fixed-point iteration in 2D.

For example, you have a Lame-Emden equation in the integral form $$h(x)=1+\int _0^R\frac{y^2}{\max (x,y)}\,\,h(y)\,y\,\text{dy}$$ where $$R$$ is a positive constant, e.g. star radius. Namely, the smallest real root of $$h(R)=0$$. The problem can be discretized using the grid $$x_i= \left(i-\frac{1}{2}\right)\text{\Delta x},\,\,\,i=1,\text{...},n$$

where $$\Delta x=\frac{R}{n}$$ is the grid spacing. Using the mid-point rule, integration amounts to matrix multiplication, i.e. our integral equation can be written as $$h_i=1+\sum _{j=1}^n M_{\text{ij}} h_j$$ where $$h_i=h (x_i)$$ and $$M_{\text{ij}}=\text{\Delta x}\text{ }\left(\frac{x_j^2}{\max \left(x_i,x_j\right)}-x_j\right)$$ Therefore, we can solve our equation by iterating the matrix as follows $$h_{\text{new}}=1+M\cdot h_{\text{old}}$$ As an initial guess, use $$h_i=h (x_i)=1$$

n=32; R=Pi; step=R/n; (*inputs*)
(* equidistant x-grid; mid-point rule allows to avoid origin, i.e. division by 0*)
x=Table[(i-0.5)*step, {i,1,n}];
(*We conversed the integral eq. into discrete summation. The scheme relies on matrix multiplication.*)
M=step*Table[(x[[j]])^2/Max[x[[i]],x[[j]]]-x[[j]], {i,1,n}, {j,1,n}];(*matrix updated at each step*)
(*we start from vector of 1's; next the vector of solutions at each grid is updated*)
h=Table[1., {i,1,n}];
err = ConstantArray[1, {100}]; (*for pre-allocation*)
count = ConstantArray[1, {100}];

eps = 10^(-14);
MaxIter = 100;
k=2;

While[k<= MaxIter && err[[k-1]] > eps, (*impose conditions*)
hold=h; (*start*)
hnext = 1+M . hold; (*dot vector onto matrix to update*)
err[[k]] = Norm[hnext-hold]; (*errors*)
h=hnext;
count[[k]]=k;
k++];
(*To display results*)
errors=Transpose[{count=DeleteCases[count,1]-1,DeleteCases[err, 1]}];

We can plot obtained numerical estimation of the root whose value approximates to the $$\pi$$ value, as it is demonstrated in the image below. Let's plot the solution
ListPlot[Transpose[{x,h}],