How to solve a non-linear system of difference equations ?
1 view (last 30 days)
Show older comments
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
2 Comments
Torsten
on 2 Jul 2019
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 ?
Answers (2)
Torsten
on 2 Jul 2019
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.
6 Comments
Torsten
on 3 Jul 2019
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.
See Also
Categories
Find more on Symbolic Math Toolbox in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!