Plotting an iterative system of nonlinear equations using MATLAB

Question

Consider the following coupled system:

Let $x(n) = \left[ \begin{array}{c} x_1(n)\\ . \\ .\\ .\\ x_{32}(n) \end{array} \right]$, and the system of $32$ first order nonlinear equations be defined by

$x_1(n + 1) = 0.2f(x_2(n)) + 0.6f(x_1(n)) + 0.2f(x_p(n)),$

$x_i(n + 1) = 0.2f(x_{i-1}(n)) + 0.6f(x_i(n)) + 0.2f(x_{i+1}(n)),$, for $2 \leq i \leq 31$, and

$x_p(n + 1) = 0.2f(x_{p-1}(n)) + 0.6f(x_p(n)) + 0.2f(x_1(n)).$

with $f(x) = 1 - \alpha x^2$ where $\alpha$ is a generic value to be experimented with.

The first part of the problem asks to find a constant matrix $A$ and column vector $g(x(n))$ so that the above system can be rewritten into a matrix-vector form as $x(n + 1) = Ag(x(n))$.

I have already found both $A$ and $x(g(n))$, they are respectively:

$A = \left[ \begin{array}{ccccccc} 0.6 & 0.2 & 0 & ..... & 0 & 0 & 0.2 \\ 0.2 & 0.6 & 0.2 & 0 &..... & 0 & 0 \\ .\\ .\\ .\\ 0 & 0 & ..... & 0 & 0.2 & 0.6 & 0.2\\ 0.2 & 0 &.....& 0 & 0 & 0.2 & 0.6 \end{array} \right]$

$g(x(n)) = \left[ \begin{array}{c} 1 - \alpha (x_1(n))^2 \\ 1 - \alpha (x_2(n))^2 \\ .\\ .\\ .\\ 1 - \alpha (x_{31}(n))^2 \\ 1 - \alpha (x_{32}(n))^2 \\ \end{array} \right] $

Now for the follow up part of this problem I am supposed to use random initial conditions and iterate the above system of difference equations and then plot $x(n)$ as a function of $n$ for three different experimental values of $\alpha = 0.5,1,$ and $1.9$. I've created an algorithm for a MATLAB code to do so, here are the steps:

1. Create constant $32\times32$ matrix A

2. Create a $32\times1$ column vector $v$ of random initial points

3. Define a array to store the different iterations of $v$

4. Create a loop iterating how the above system is defined so $x(n + 1) = Ag(x(n))$

Well I have managed to succeed in coding steps 1-3 and some of step 4, here is the code I have so far:

A = [0.6 0.2];

B = zeros(1,29);

C = [0.2];

C = [A B C];

i = 1;

D = circshift(C,1);

E = [C; D];

while i < 16

C = circshift(D,1);

E = [E; C];

D = circshift(C,1);

E = [E; D];

i = i + 1;

end

alpha = 0.5;

x= zeros(32,1);

v = rand(32,1);

for n = 2:1000

x(n) = E*(1 - alpha.*v(n - 1).^2);

end

plot(1:2000,x);

I'm trying to finish the rest of the code but need help, the final loop I have in the above code for completing step 4 of my algorithm is incomplete but the rest of the code should work fine. Any help with this would be appreciated, thanks.

Answer 1

The way I see it, you will need to define the function $g$ that takes as argument a vector and returns another vector, something like:

function y = fun_g(x, alpha)

y = 1 - alpha .* x .* x; % y will be a vector of the same size as x

And in your main function add near the end something like:

nmax = 100;
x_t = zeros(32,nmax);
n = 1;
x_t(:,1) = x;

while (n <= nmax)
    x_t(:,n+1) = E * fun_g(x_t(:,n), alpha);
    n = n+1;
end

And then do the plots.

I'm sorry I haven't tested my code, so I might have missed something.

---

By the way, I assumed that matrix $A$ in your post is matrix $E$ in your code. I'm not sure because the way you created matrix $E$ was confusing for me.

There are easier ways of creating a toeplitz matrix, with the Matlab reference.