Solving non linear system with fsolve

I'm trying to solve a system of non linear equations:

Each element of the matrix R is multiplied by the boltzmann constant so it is extremely small, while C contains elements that are much bigger. When I try to solve the system with fsolve i receive the following error and the solution makes no physical sense.

This is my code:

R = 1.0e-20 * [
    0.1251   -0.0001   -0.0000   -0.1019   -0.0008
   -0.0001    0.4489   -0.0223   -0.0210   -0.4025
   -0.0000   -0.0223    0.1251   -0.0008   -0.1019
    0.1019    0.0210    0.0008   -0.6114    0.0852
    0.0008    0.4025    0.1019    0.0852   -0.6114];
C = 1.0e+05 * [

    0.3600         0         0   -0.3600         0
         0    1.4400         0         0   -1.4400
         0         0    0.3600         0   -0.3600
    0.3600    1.4400         0   -2.2958         0
         0    1.4400    0.3600   -0.4958   -2.2958];
f=[
1.67220567523920e-11
2.14091795584358e-12
9.87213173364192e-14
-152843740.074950
-1.57644730491682e-11];

Function= @(Temp) R*Temp.^4+C*Temp-f;
TGuess=800*ones(5,1);
Tsolution=fsolve(Function,Tguess);

Does anybody have suggestions on how to avoid the stalling of the solver?

[Walter Roberson]

You have multinomials of degree 4. Each variable has 4 solution families, and there are 5 variables, so you have 4^5 = 1024 solutions. Some of the solutions might involve complex components.

Answer 1

[Star Strider]

It seems that fsolve stopped when it converged.

R = 1.0e-20 * [
    0.1251   -0.0001   -0.0000   -0.1019   -0.0008
   -0.0001    0.4489   -0.0223   -0.0210   -0.4025
   -0.0000   -0.0223    0.1251   -0.0008   -0.1019
    0.1019    0.0210    0.0008   -0.6114    0.0852
    0.0008    0.4025    0.1019    0.0852   -0.6114];
C = 1.0e+05 * [
    0.3600         0         0   -0.3600         0
         0    1.4400         0         0   -1.4400
         0         0    0.3600         0   -0.3600
    0.3600    1.4400         0   -2.2958         0
         0    1.4400    0.3600   -0.4958   -2.2958];
f=[
1.67220567523920e-11
2.14091795584358e-12
9.87213173364192e-14
-152843740.074950
-1.57644730491682e-11];
Function= @(Temp) R*Temp.^4+C*Temp-f;
TGuess=800*rand(5,1);
format longG
[Tsolution,fcnval]=fsolve(Function,TGuess)

Equation solved, solver stalled.

fsolve stopped because the relative size of the current step is less than the
value of the step size tolerance squared and the vector of function values
is near zero as measured by the value of the function tolerance.

Tsolution = 5×1

452.763019358226 -452.763019358226 -452.763019358226 452.763019358226 -452.763019358226

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

fcnval = 5×1

1.0e+00 * -7.35099816949493e-12 -8.80237428996419e-13 -5.66986331081861e-14 2.98023223876953e-08 6.9397093612381e-12

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

The function values at convergence are appropriately very small.

.

Answer 2

[Matt J]

R is too small compared to C to have any effect numerically, so the equation is basically C*T=f. I'm assuming also that you meant to constrain it so that the temperature are positive, T>=0. That means you should use lsqnonneg instead,

R = 1.0e-20 * [
    0.1251   -0.0001   -0.0000   -0.1019   -0.0008
   -0.0001    0.4489   -0.0223   -0.0210   -0.4025
   -0.0000   -0.0223    0.1251   -0.0008   -0.1019
    0.1019    0.0210    0.0008   -0.6114    0.0852
    0.0008    0.4025    0.1019    0.0852   -0.6114];

C = 1.0e+05 * [
    0.3600         0         0   -0.3600         0
         0    1.4400         0         0   -1.4400
         0         0    0.3600         0   -0.3600
    0.3600    1.4400         0   -2.2958         0
         0    1.4400    0.3600   -0.4958   -2.2958];
f=[
1.67220567523920e-11
2.14091795584358e-12
9.87213173364192e-14
-152843740.074950
-1.57644730491682e-11];

R=R*1e-5; C=C*1e-5; f=f*1e-5; %changes nothing in infinite-precision math

Temp0=lsqnonneg(C, f)

Temp0 = 5×1

703.2679 75.3289 405.0601 807.0125 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

We can examine what would happen if we included the R term using lsqnonlin. Essentially nothing, as it turns out:

Func= @(Temp) 1e25*R*Temp.^4+1e25*C*Temp-1e25*f;

Temp=lsqnonlin(Func,800*ones(5,1), zeros(5,1))

Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

Temp = 5×1

703.2680 75.3289 405.0601 807.0126 0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

norm(Temp-Temp0,inf)

ans = 6.8752e-05

Answer 3

[John D'Errico]

Whenever you have a problem with hugely varying constants in it like this, it makes sense to consider if you should transform things in a subtle way.

For example, here, I will factor out the 1e-20 from R. Don't worry. I'll remember it was there. But to make this work in double precision arithmetic, we will need this. I could also pull the 1e5 from C. Don't worry. (Be happy.)

R = [0.1251   -0.0001   -0.0000   -0.1019   -0.0008
    -0.0001    0.4489   -0.0223   -0.0210   -0.4025
    -0.0000   -0.0223    0.1251   -0.0008   -0.1019
    0.1019    0.0210    0.0008   -0.6114    0.0852
    0.0008    0.4025    0.1019    0.0852   -0.6114];
C = [0.3600         0         0   -0.3600         0
    0    1.4400         0         0   -1.4400
    0         0    0.3600         0   -0.3600
    0.3600    1.4400         0   -2.2958         0
    0    1.4400    0.3600   -0.4958   -2.2958];

T will be the vector of unknowns. transform the problem like this:

R*(T*1e-5)^4 + C*(T*1e5) = f

Do you see the problem is identical to yours? We are internally raising 1e-5 to the 4th power. And that will make the numerical issues simpler. I have put those nasty powers of 10 in a different place. Now the problem MAY be more easily solved using double precision. EXCEPT...

A numerical problem will still arise though. I think you messed up on one of the elements of f.

f=[ 1.67220567523920e-11
    2.14091795584358e-12
    9.87213173364192e-14
    -152843740.074950
    -1.57644730491682e-11];

Do you see that f(4) is insanely different in magnitude from the rest? It differes by 20 powers of 10. And again, double precision arithmetic will fail miserably. I think it is unlikely to be true that is the correct value of f(4).

Until you fix that problem, and tell us what number really should live at f(4), you can go no further.

[Riccardo]

I'm double checking the passages that lead to those matrices, but i think that the difference in the order of magnitude in f should be correct. I'm trying to find the temperatures of gasses and walls in a thermal plant combustor. Since the equations are related to the energy balances of the system and since the fourth line is related to the combusiton chamber where there is the production of 70 MW of energy, it makes sense that f(4) is much bigger than the other terms