How to solve a non-linear system of difference equations ?
4 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Torsten
le 2 Juil 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 ?
Réponses (2)
Torsten
le 2 Juil 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 commentaires
Torsten
le 3 Juil 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.
Voir également
Catégories
En savoir plus sur Symbolic Math Toolbox dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!