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?

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

3 Answers 3

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

enter image description here

  • Calculate the fixed point by DiscrteLimit.
 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 -> ∞]




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;

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*)
(*To display results*)
errors=Transpose[{count=DeleteCases[count,1]-1,DeleteCases[err, 1]}];
head = {"Iterations", "Error"};
Print["The solution vector at each point of grid is: ", h]

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

         PlotStyle->Magenta, AxesLabel->{"x","h(x)"}, PlotLabel->"Numerical solution" ,

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.