How to solve a non-linear system of difference equations ?

Hi ! I am a beginner with matlab and I have some promblem very important for my research...

I am trying to solve (or more precisely to simulate) the following system numerically.

I write the following code but I think it misses some matlab function (like solve for exemple)...

The point it that I don't know how I can implement it together with the loop...

Here is the code :

K(1) = 10;
S(1) = 1000;
S(2) = 950;
E(1) = 1000;
A(1) = 10;
a = 0.9;
alpha = 0.3;
beta = 0.6;
gamma = 0.1;
z = 0.2;
rho = 2;
delta = 0.9;
phi = 0.2;
b = 0.1;
Etilde = 1000;
%b*Etilde+(1-b)*E(t);

x = linspace(1,100);
for t=1:length(x)-2
    S(t+1) = (1/(1+rho))*(((1+alpha*A(t+1)*K(t+1)^(alpha-1)*(1+z)^beta*(S(t+1)-S(t+2))^gamma-delta)*beta*A(t)*K(t)^alpha*(1+z)^beta*(S(t)-S(t+1))^gamma)/((1+alpha*A(t+1)*K(t+1)^(alpha-1)*(1+z)^beta*(S(t+1)-S(t+  2))^gamma-delta)*gamma*A(t)*K(t)^alpha*(1+z)^beta*(S(t)-S(t+1))^(gamma-1)-gamma*A(t+1)*K(t+1)^alpha*(1+z)^beta*(S(t+1)-S(t+2))^(gamma-1)));
    K(t+1) = beta*A(t)*K(t)^alpha*(1+z)^beta*(S(t)-S(t+1))^gamma-gamma*A(t+1)*K(t+1)^alpha*(1+z)^beta*(S(t+1)-S(t+2))^(gamma-1)*S(t+1);
    E(t+1) = (1/(1+rho))*(b*Etilde+(1-b)*E(t)+phi*((((1+alpha*A(t+1)*K(t+1)^(alpha-1)*(1+z)^beta*(S(t+1)-S(t+2))^gamma-delta)*beta*A(t)*K(t)^alpha*(1+z)^beta*(S(t)-S(t+1))^gamma)/((1+alpha*A(t+1)*K(t+1)^(alpha-1)*(1+z)^beta*(S(t+1)-S(t+2))^gamma-delta)*gamma*A(t)*K(t)^alpha*(1+z)^beta*(S(t)-S(t+1))^(gamma-1)-gamma*A(t+1)*K(t+1)^alpha*(1+z)^beta*(S(t+1)-S(t+2))^(gamma-1)))-S(t)));
    A(t+1) = A(t)*(1+a);
end
K(t)
E(t)
S(t)
figure(1)
plot(x,K)
figure(2)
plot(x,E)
figure(3)
plot(x,S)
grid

If someone has an idea...

Many thanks,

Nicolas

[Torsten]

Most probably, you won't be able to get an analytical formula for S, T and E depending on the initial values and the index t. Or what did you expect to get ?

[Nicloot]

Sure, that's why I would like to simulate the model in order to plot K, S and E for t=0...N (N large...). The point is that it would require some numerical solver, but I don't know which one I should use (vpasole???) and how to do it.

Answer 1

[Torsten]

Call "fsolve" repeatedly to get

(S3,K2,E2),
(S4,K3,E3),
...
(S(t+2),K(t+1),E(t+1))

A(t) = A(1)*(1+a)^(t-1) 

is obvious.

[Nicloot]

Actually, this is the point : I tried to use fsolve (and other solvers as vpasolve) inside the loop but I think this is not the good way to do. As I wrote in the question, I am really not familiar with matlab.

Anyway, thank you for your answer, I will try to iterate fsolve !

[Torsten]

K = 10;
So = 1000;
S = 950;
E = 1000;
A = 10;
a = 0.9;
alpha = 0.3;
beta = 0.6;
gamma = 0.1;
z = 0.2;
rho = 2;
delta = 0.9;
phi = 0.2;
b = 0.1;
Etilde = 1000;
for t=1:100
  x0 = [S;K;E;A];
  fun = @(x)[S-( (1/(1+rho))*(((1+alpha*x(4)*x(2)^(alpha-1)*(1+z)^beta*(S-x(1))^gamma-delta)*beta*A*K^alpha*(1+z)^beta*(So-S)^gamma)/((1+alpha*x(4)*x(2)^(alpha-1)*(1+z)^beta*(S-x(1))^gamma-delta)*gamma*A*K^alpha*(1+z)^beta*(So-S)^(gamma-1)-gamma*x(4)*x(2)^alpha*(1+z)^beta*(S-x(1))^(gamma-1)));
x(2) - (beta*A*K^alpha*(1+z)^beta*(So-S)^gamma-gamma*x(4)*x(2)^alpha*(1+z)^beta*(S-x(1))^(gamma-1)*S);
x(3) - ( (1/(1+rho))*(b*Etilde+(1-b)*E+phi*((((1+alpha*x(4)*x(2)^(alpha-1)*(1+z)^beta*(S-x(1))^gamma-delta)*beta*A*K^alpha*(1+z)^beta*(So-S)^gamma)/((1+alpha*x(4)*x(2)^(alpha-1)*(1+z)^beta*(S-x(1))^gamma-delta)*gamma*A*K^alpha*(1+z)^beta*(So-S)^(gamma-1)-gamma*x(4)*x(2)^alpha*(1+z)^beta*(S-x(1))^(gamma-1)))-So)));
x(4)-A*(1+a)];
  sol = fsolve(fun,x0)
  So = S;
  S = sol(1);
  K = sol(2);
  E = sol(3);
  A = sol(4);
end 

I leave it to you to save the complete solution.

[Nicloot]

Many thanks ! I wrote the following, still a problem...

clc;clear;

K = 1;
So = 1000;
S = 950;
E = 1000;
A = 1;
a = 0.91;
alpha = 0.3;
beta = 0.6;
gamma = 0.1;
z = 0.2;
rho = 2;
delta = 0.9;
phi = 0.2;
b = 0.1;
Etilde = 10000;

S_opt=NaN(1,100);
K_opt=NaN(1,100);
E_opt=NaN(1,100);
A_opt=NaN(1,100);


for t=1:100
  x0 = [S;K;E;A];
  fun = @(x)[S-( (1/(1+rho))*(((1+alpha*x(4)*x(2)^(alpha-1)*(1+z)^beta*(S-x(1))^gamma-delta)...
      *beta*A*K^alpha*(1+z)^beta*(So-S)^gamma)/((1+alpha*x(4)*x(2)^(alpha-1)*(1+z)^beta*(S-x(1))^gamma-delta)...
      *gamma*A*K^alpha*(1+z)^beta*(So-S)^(gamma-1)-gamma*x(4)*x(2)^alpha*(1+z)^beta*(S-x(1))^(gamma-1)))...
      -((rho)/(1+rho))*((b*Etilde+(1-b)*E)/(phi)-So));
      
    x(2) - (beta*A*K^alpha*(1+z)^beta*(So-S)^gamma-gamma*x(4)*x(2)^alpha*(1+z)^beta*(S-x(1))^(gamma-1)*S);
    
    x(3) - ( (1/(1+rho))*(b*Etilde+(1-b)*E+phi*((((1+alpha*x(4)*x(2)^(alpha-1)...
        *(1+z)^beta*(S-x(1))^gamma-delta)*beta*A*K^alpha*(1+z)^beta*(So-S)^gamma)/...
        ((1+alpha*x(4)*x(2)^(alpha-1)*(1+z)^beta*(S-x(1))^gamma-delta)*gamma*A*K^alpha...
        *(1+z)^beta*(So-S)^(gamma-1)-gamma*x(4)*x(2)^alpha*(1+z)^beta*(S-x(1))^(gamma-1)))-So)));
        
    x(4)-A*(1+a)];

  sol = fsolve(fun,x0);
  So = S;
  S = sol(1);
  K = sol(2);
  E = sol(3);
  A = sol(4);
  
  S_opt(t)= sol(1);
  K_opt(t) = sol(2);
  E_opt(t) = sol(3);
  A_opt(t) = sol(4);
end 


T=1:1:100;

figure(1)
plot(T,K_opt)
figure(2)
plot(T,E_opt)
figure(3)
plot(T,S_opt)
grid

matlab answers :

Error using trustnleqn (line 28)

Objective function is returning undefined values at initial point. FSOLVE cannot continue.

Error in fsolve (line 408)

trustnleqn(funfcn,x,verbosity,gradflag,options,defaultopt,f,JAC,...

Error in myprogr(line 41)

sol = fsolve(fun,x0);

[Torsten]

For which t value does the message appear ?

[Nicloot]

Immediately in the first loop. I try several combinations of parameter values and initial conditions, but they all fail...

[Torsten]

You can test if "fun" returns reasonable values if you insert the line

fun(x0)

before calling "fsolve".

I think you will come into trouble with your exponentiations x^a if x becomes negative or if x becomes 0 and a is negative.

Answer 2

[Nicloot]

Here is the matlab output,

ans =

1.0e+03 *

3.6167

Inf

0.7333

-0.0009

There is clearly a problem since we can see that x(4)=-0.9 which is clearly different than the expected 1.91 (since)

(edit : I change the parameter a to 0.91 and initial value for A to 1 in my last attempt)

[Torsten]

The output gives the function value, not the variable value:

x(4)-x0(4)*(1+a) = -0.9

The bigger problem is fun(2) = Inf which results from 0^(gamma-1)

[Nicloot]

Ok, thanks, so the value for A is logical.

I checked my calculation twice... Don't understand what goes wrong... I'll check one more time !

Thank again !