How to use fsolve within multiple for loops?

10 次查看（过去 30 天）
Ronel Scheepers
Ronel Scheepers 2020-1-13
评论： Ronel Scheepers 2020-1-13
I have a system of 4 non-linear equations with 4 unknowns x(1), x(2), x(3) and x(4). I am able to solve the equations using the following code, however, fsolve only gives me one possible solution. I determined analytically that there is one positive steady state but would like find the approximation to the SS numerically using loops for the initial guesses as well as for different parameter values and are struggling with the code for the for loops. Any help appreciated.
% The function accepts as input the variables
%(x, m0 ,m1, gamma, C0, r1, r2, Uc, p, r4, r5, r6, a, b, r7),
% where x = [x(1) x(2) x(3) x(4)] = [C U E M],
% and returns vector F(x) as output.
%Parameters taken as input
m0=0.1; m1=1; gamma=1; C0=1; r1=0.01;r2=0.01; Uc=0.25; p=1.0; r4=1.0;
r5=0.01; r6=1.0;a=0.1; b=0.1;r7=0.1;
% Calling fsolve
fun = @(x)SteadyStates(x,m0 ,m1, gamma, C0, r1, r2, Uc, p, r4, r5, r6, a, b, r7);
x0 = [0.5;0.5;0.2;0.2]; % initial values
x = fsolve(fun,x0)
function F = SteadyStates(x, m0 ,m1, gamma, C0, r1, r2, Uc, p, r4, r5, r6, a, b, r7)
C=x(1);
U=x(2);
E=x(3);
M=x(4);
F(1) = m0*(1+m1*gamma*x(1))/(1+gamma*x(1))*(C0-x(1))-r1*x(1);
F(2) = r1*x(1)+r2*x(2)*(Uc-x(2))-x(2)^2/(p+x(2))+r4*x(3)/(1+x(2))-r5*x(2)-r6*x(2)*x(3);
F(3) = a*r1*x(1)+x(2)^2/(p+x(2))-r4*x(3)/(1+x(2))-b*r6*x(2)*x(3);
F(4) = b*r6*x(2)*x(3) - r7*x(4);
end

请先登录，再回答此问题。

采纳的回答

Walter Roberson
Walter Roberson 2020-1-13
If you have the Symbolic Toolbox, then solve F(1), F(2), F(4) for x(1), x(3), x(4) . substitute those values into F(3) and simplify, to get
(4*x2^4 + 4407*x2^3 + 43*x2^2 + 30*x2 - 40)/(4000*(x2^2 + x2 - 1))
so three of the variables can be be written in terms of x(2) which in turn gets solved as a quartic (degree 4)
You can get the complete detailed solution by solve(F,'MaxDegree',4) but I find those kinds of solutions to be essentially useless.
You can also use vpasolve() on the set of four equations. If you do that then the result is 5 solutions and a warning about possible numeric accuracy problems with one of them. But the first of the 5 solutions is a false solution that involves division by 0. You can demonstrate that the theoretical solution corresponding to the one it complains about accuracy for, is in fact an exact solution.
  1 个评论
Ronel Scheepers
Ronel Scheepers 2020-1-13
Thank you so much, vpasolve() worked. I now have to see if I can include a loop to vary the input of parameter values as the current values is just a randomly chosen set.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by