Trouble getting steady-state solution by solving system of nonlinear algebraic equations in MATLAB

Question

Background

I have a stiff system of 6 ODEs, represented in MATLAB as follows:

system = @(t,x)[x(1,:).*x(2,:).*(-5.726882618492327e8)-x(1,:).*x(3,:).*1.467710449114545e10-x(1,:).*x(4,:).*3.012162507288153e12+x(3,:).*x(6,:).*1.137674873581712e12+x(4,:).*x(6,:).*5.703626484852603e12+x(5,:).*x(6,:).*1.656793950858653e9-x(1,:).^2.*2.980001443206009e8+x(6,:).^2.*8.25e9;x(1,:).*x(2,:).*(-5.726882618492327e8)+x(3,:).*x(6,:).*1.137674873581712e12-x(2,:).^2.*4.476510081133096e2+x(6,:).^2.*3.403491297089207e3;x(1,:).*x(2,:).*5.726882618492327e8-x(1,:).*x(3,:).*1.467710449114545e10-x(3,:).*x(6,:).*1.137674873581712e12+x(4,:).*x(6,:).*5.703626484852603e12;x(1,:).*x(3,:).*1.467710449114545e10-x(1,:).*x(4,:).*3.012162507288153e12-x(4,:).*x(6,:).*5.703626484852603e12+x(5,:).*x(6,:).*1.656793950858653e9;x(5,:).*(-8.186979330234089e8)+x(6,:).*2.1e9+x(1,:).*x(4,:).*3.012162507288153e12-x(5,:).*x(6,:).*1.656793950858653e9;x(5,:).*8.186979330234089e8-x(6,:).*2.1e9+x(1,:).*x(2,:).*5.726882618492327e8+x(1,:).*x(3,:).*1.467710449114545e10+x(1,:).*x(4,:).*3.012162507288153e12-x(3,:).*x(6,:).*1.137674873581712e12-x(4,:).*x(6,:).*5.703626484852603e12-x(5,:).*x(6,:).*1.656793950858653e9+x(1,:).^2.*2.980001443206009e8+x(2,:).^2.*4.476510081133096e2-x(6,:).^2.*8.250003403491297e9];

If I use an ODE solver, I can solve the ODE (shown here) and use it to approximate the steady-state solution. As an example,

[t,x] = ode15s(system,[0 1e-6],[0 0 0 0 0 1]);
plot(t,x)
ss_sol = x(end,:);

For reference, the steady-state solution obtained in this way is approximately

ss_sol = [0.322330352943109, 0.458043435766086, 0.001213186698607, 0.000016426443105, 0.157136142002196, 0.061260531336451];

Even if I change the initial conditions a bit, I come back to this solution. Note that sum(ss_sol) is approximately 1 (the ODEs were derived under this assumption). Also note that due to physical constraints, each variable must be between 0 and 1, which the solution satisfies.

Problem

To get the steady-state solution, I'd like to set the derivatives equal to zero and solve the analogous, nonlinear algebraic system for its roots. The system is slightly rewritten in MATLAB as

system = @(x)[x(1).*x(2).*(-5.726882618492327e8)-x(1).*x(3).*1.467710449114545e10-x(1).*x(4).*3.012162507288153e12+x(3).*x(6).*1.137674873581712e12+x(4).*x(6).*5.703626484852603e12+x(5).*x(6).*1.656793950858653e9-x(1).^2.*2.980001443206009e8+x(6).^2.*8.25e9,x(1).*x(2).*(-5.726882618492327e8)+x(3).*x(6).*1.137674873581712e12-x(2).^2.*4.476510081133096e2+x(6).^2.*3.403491297089207e3,x(1).*x(2).*5.726882618492327e8-x(1).*x(3).*1.467710449114545e10-x(3).*x(6).*1.137674873581712e12+x(4).*x(6).*5.703626484852603e12,x(1).*x(3).*1.467710449114545e10-x(1).*x(4).*3.012162507288153e12-x(4).*x(6).*5.703626484852603e12+x(5).*x(6).*1.656793950858653e9,x(5).*(-8.186979330234089e8)+x(6).*2.1e9+x(1).*x(4).*3.012162507288153e12-x(5).*x(6).*1.656793950858653e9,x(5).*8.186979330234089e8-x(6).*2.1e9+x(1).*x(2).*5.726882618492327e8+x(1).*x(3).*1.467710449114545e10+x(1).*x(4).*3.012162507288153e12-x(3).*x(6).*1.137674873581712e12-x(4).*x(6).*5.703626484852603e12-x(5).*x(6).*1.656793950858653e9+x(1).^2.*2.980001443206009e8+x(2).^2.*4.476510081133096e2-x(6).^2.*8.250003403491297e9];

However, I have trouble getting the same values of ss_sol mentioned earlier.

Using vpasolve, even when I increase the precision to 100+ digits, does not lead to any solution. Using fsolve or lsqnonlin (the latter of which lets me set the lower and upper bounds of 0 and 1, respectively) converge to values that are not the steady-state solution, unless I supply ss_sol as the initial guess. I try decreasing the tolerances and increasing the number of iterations but still no luck.

Any suggestions? Is the best way really to just integrate the ODEs to steady-state? I keep reading that reducing the stiff system of ODEs to an algebraic equation should be easier to solve, but I just can't find the steady-state solutions this way. I need a general way to solve for the steady-state values for many different systems, so the problem extends beyond the aforementioned specific example.

Answer 1

The problems you are encountering are not unusual.

Nonlinear systems of equations, as you have already discovered, often have multiple solutions. It can be difficult to find the specific one you want.

Basic Newton's method is not guaranteed to converge to any solution given a poor initial guess. (I had thought that the MATLAB/fsolve algorithm would be able to find some solution in this case, however.)

There is a huge body of literature on techniques for solving nonlinear algebraic equations. A large subset of this literature deals with algorithms called "continuation methods" designed to deal with particularly challenging problems. Many of these continuation methods are very similar to ODE solution methods!

At this point, I'll assume you are more interested in simply solving your equations than researching numerical methods.

Here is what I suggest:

If you have a good initial guess, basic Newton's methods will converge. In your case, it looks like your ode15s solution is approximately steady-state at t=2e-7 which is significantly before your actual termination time of t=1e-6. Since you only need an approximate guess, you can probably also reduce the accuracy of ode15s (by setting the options AbsTol and RelTol) by at least an order of magnitude. Try calculating an initial guess by this approach and then passing it to fsolve to calculate the final steady-state solution.