A Solve Badly Scaled Problem with Newton-Chord Technique

[Ulver48]

Hello guys,

I try to use the Newton - Chord technique in order to solve a nonlinear system and find it's equilibrium points.This method requires the inverse of the Jacobian Matrix of the nonlinear system. After the linearization around the given starting point x₀, I create a linear approximation of the initial problem and then I compute the Jacobian Matrix. Let's call this matrix J. The problem is that this matrix is badly scaled. The system of ODEs consists of 5 equations, 3 of which describe the rate of change for some probabilities P1,P2,P3 (so the range of their values are somewhere between 0 and 1) and two equations that describe the rate of change for an Electric Field ( the real and the imaginary part, who are normalized in such a way that they take values of the order of 1e17).

I use Matlab which tells me that the resulting Jacobian is very badly scaled ( for instance rcond = 1e-18). I tried replacing the Electric Field by a scaled version F, which is equal to the electric field divided by the value 1e18, but it didn't work. So I wanted to ask for your advice. What would someone do in a situation like this ? If you want further information then I will give you more details about my particular set of equations and the values of the resulting Jacobian. I wanted to avoid a really long post and I know that most of you are busy with your own work. Time is a precious thing.I hope for some general guidelines.that could help me deal with my problem. Thank you all for your time.

 

[jedishrfu]

I think you're going to have to sketch out what you've done so some of our mentors can look at it and comment. From what you've written its too fuzzy to really give advice. Also it might require a mentor to type in the equations that you didn't in order to explain what should be done.

Having said that, you should also be aware that you may not get a timely answer as we are all volunteers here and tackle thread questions when we can.

Not understanding the problem you are trying to solve, I have to ask by scaling you mean that some values of the matrix are very very large and/or some values are very very small, right?

Matlab has a website where students ask these kinds of questions often and I found these responses:

https://www.mathworks.com/matlabcen...atrix-is-singular-or-badly-scaled-please-help

https://www.mathworks.com/matlabcen...g-matrix-is-close-to-singular-or-badly-scaled

and there are other pages I didn't explore that you can:

https://www.google.com/search?q=matlab+how+to+handle+badly+scaled+matrices

 

[Ulver48]

Thank you very much for your response. My set of equations is
$$
\dfrac{dP_{wl}}{dt}= 5.8858*10^8+2.2071*10^8*P_{es}-4.0833*10^{9}*P_{wl}+1.8625*10^9*P_{es}*P_{wl} \\
\dfrac{dP_{gs}}{dt}= -1.9514*10^{-10}*(E_{real}^2+E_{imag}^2)*(2*P_{gs}-1)+3.75*10^{11}*P_{es}*(1-P_{gs})+5.5556*10^{10}*P_{gs}*(1-P_{es})-8.3333*10^8*P_{gs}\\
\dfrac{dP_{es}}{dt} = 9.439*10^{11}*P_{wl}*(1-P_{es})+1.6667*10^{11}*P_{gs}*(1-P_{es})-1*10^{11}*P_{es}*(1-P_{wl})-1.25*10^{11}*P_{es}*(1-P_{gs})\\
\dfrac{dE_{real}}{dt} = 0.5*(6.0979*10^{10}*(2Pgs-1)-6.3104*10^{10})*(E_{real}-1.3*E_{imag})-\Delta\omega*E_{imag}+k_{eff}*1.8269*10^{19}\\
\dfrac{dE_{imag}}{dt} = 0.5*(6.0979*10^{10}*(2Pgs-1)-6.3104*10^{10})*(E_{imag}+1.3*E_{real})+\Delta\omega*E_{real}
$$

The terms $k_{eff},\Delta\omega$ are terms that I can control. The Δω term is equalt to 2πf , where f is the frequency (somewhere around GHz). Now I compute the Jacobian with respect to some initial values for my co-ordinates. The Jacobian is quite a monster, so I will post my Matlab code.

%%%%%%%%%%%%%%%%%%HERE COMES THE CODE%%%%%%%%%%%%%%%%%%%%%
% Fwl - Terms
dFwldPwl = -4.0833e9+Pes*1.8625e9;
dFwldPgs=0;
dFwldPes=2.2071e8+Pwl*1.8625e9;
dFwldEreal = 0;
dFwldEimag=0;

%Fgs - Terms
dFgsdPwl=0;
dFgsdPgs=-1.9514e-10*2*(Ereal^2+Eimag^2)-3.75e11*Pes-8.3333e8-5.5556e10*(1-Pes);
dFgsdPes=3.75e11*(1-Pgs)+5.5556e10*Pgs;
dFgsdEreal=-1.9514e-10*(2*Ereal)*(2*Pgs-1);
dFgsdEimag=-1.9514e-10*(2*Eimag)*(2*Pgs-1);

%Fes - Terms
dFesdPwl=9.439e11*(1-Pes)+1e11*Pes;
dFesdPgs=1.6667e11*(1-Pes)+1.25e11*Pes;
dFesdPes=-9.439e11*Pwl-1.6667e11*Pgs-1e11*(1-Pwl)-1.25e11*(1-Pgs)-2e9;
dFesdEreal=0;
dFesdEimag=0;

%F_Ereal - Terms
dF_ErealdPwl=0;
dF_ErealdPgs=0.5*6.097e10*2*(Ereal-1.3*Eimag);
dF_ErealdPes=0;
dF_ErealdEreal=0.5*(6.0979e10*(2*Pgs-1)-6.3104e10);
dF_ErealdEimag=0.5*(6.0979e10*(2*Pgs-1)-6.3104e10)*(-1.3)-Delta;

%F_Eimag - Terms
dF_EimagdPwl=0;
dF_EimagdPgs=0.5*6.097e10*2*(Eimag+1.3*Ereal);
dF_EimagdPes=0;
dF_EimagdEreal=0.5*(6.0979e10*(2*Pgs-1)-6.3104e10)*(1.3)+Delta;
dF_EimagdEimag=0.5*(6.0979e10*(2*Pgs-1)-6.3104e10);

J=[dFwldPwl dFwldPgs dFwldPes dFwldEreal dFwldEimag;
    dFgsdPwl dFgsdPgs dFgsdPes dFgsdEreal dFgsdEimag;
    dFesdPwl dFesdPgs dFesdPes  dFesdEreal dFesdEimag;
    dF_ErealdPwl dF_ErealdPgs dF_ErealdPes dF_ErealdEreal dF_ErealdEimag;
    dF_EimagdPwl dF_EimagdPgs dF_EimagdPes dF_EimagdEreal dF_EimagdEimag];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Now for values
$$P_{wl}=0.2712, P_{gs}=0.5,P_{es}=0.7148,E_{real}=6.8158*10^{17},E_{imag}=-4.7887*10^{17},k_{eff}=1.7582*10^9,\Delta\omega=Delta=1.8850*10^{10} $$
The Jacobian is
$$ J = \begin{bmatrix}
-2.7520*10^9 & 0 &7.2583*10^8 & 0 & 0 \\
0 & -2.7080*10^{26} & 2.1528*10^{11} & -5.3159*10^{-7} & 3.7349*10^{-7} \\
3.4067*10^{11}&1.3688*10^{11} &-4.7671*10^{11} & 0 &0 \\
0& 7.9511*10^{28} &0 &-3.1552*10^{10} &2.2168*10^{10} \\
0& 2.4826*10^{28} &0 &-2.2168*10^{10} &-3.1552*10^{10} \end{bmatrix}
$$

This matrix is clearly badly scaled. The condition number is 2.2017e+24, which is very big and I have to invert this matrix.

 

[mfb]

The matrix in WolframAlpha (only works sometimes, refresh to try again if necessary). It gets a different, but still very large condition number.

I tried replacing the Electric Field by a scaled version F, which is equal to the electric field divided by the value 1e18, but it didn't work.

How does the matrix look like after this scaling?

 

[Ulver48]

First of all let me tell you how I scaled my matrix. Let's write the equations this way
$$ \dfrac{dP_{wl}}{dt}=F_{wl}(P_{wl},P_{gs},P_{es},E_{real},E_{imag})\\
\dfrac{dP_{gs}}{dt}=F_{gs}(P_{wl},P_{gs},P_{es},E_{real},E_{imag})\\
\dfrac{dP_{es}}{dt}=F_{es}(P_{wl},P_{gs},P_{es},E_{real},E_{imag})\\
\dfrac{dE_{real}}{dt}=F_{real}(P_{wl},P_{gs},P_{es},E_{real},E_{imag})\\
\dfrac{dE_{imag}}{dt}=F_{imag}(P_{wl},P_{gs},P_{es},E_{real},E_{imag})$$

Then I use two new variables $E^n_{real}=10^{-18}*E_{real}$ and $E^n_{imag}=10^{-18}*E_{imag}$. Then I substitute them in my equations above and I end up with the following form.
$$ \dfrac{dP_{wl}}{dt}=F_{wl}(P_{wl},P_{gs},P_{es},10^{18}*E^n_{real},10^{18}*E^n_{imag})\\
\dfrac{dP_{gs}}{dt}=F_{gs}(P_{wl},P_{gs},P_{es},10^{18}*E^n_{real},10^{18}*E^n_{imag})\\
\dfrac{dP_{es}}{dt}=F_{es}(P_{wl},P_{gs},P_{es},10^{18}*E^n_{real},10^{18}*E^n_{imag})\\
\dfrac{dE^n_{real}}{dt}=\dfrac{1}{10^{18}}F_{real}(P_{wl},P_{gs},P_{es},10^{18}*E^n_{real},10^{18}*E^n_{imag})\\
\dfrac{dE^n_{imag}}{dt}=\dfrac{1}{10^{18}}F_{imag}(P_{wl},P_{gs},P_{es},10^{18}*E^n_{real},10^{18}*E^n_{imag})$$

So I end up with the following Jacobian matrix.
$$J=\begin{bmatrix}
-2.7520*10^9 & 0 & 7.2583*10^8 & 0 & 0 \\
0 & -2.7080*10^26& 2.152*10^11 & -5.3159*10^{-7}& 3.7349*10^{-7}\\
3.4067*10^{11} & 1.3688*10^{11} & -4.7671*10^{11} & 0 & 0 \\
0 & 7.9511*10^{10} & 0 & -3.1552*10^{-8} & 2.2168*10^{-8}\\
0 & 2.4826*10^{10}& 0 & -2.2168*10^{-8} & -3.1552*10^{-8}
\end{bmatrix} $$

According to Matalb, the condition number of this matrix is 7.0226e+33, which is worse than the previous one. Of course I used the same values for my co-ordinates, except for the field values which I divided by $10^{18}$. To make things worse, the Newton method itself adds a vector in my matrix ( one more row at the end), which is normalized ( values between 0 and 1) before inverting it. I use a tool called matcont in order to create a bifurcation diagram with respect to the two control parameters and it uses the Newton method to correct my equilibrium point. If you don't know what a bifurcation diagram is , it's ok. My goal now is just to create a well-behaved Jacobian matrix.

 

[mfb]

Where did you take the factor 10¹⁸ into account for the last two columns? You applied it to the rows, but the 10¹⁸ in the argument has to be considered as well. In WolframAlpha that improves the condition number by 5.5 orders of magnitude.

In the derivative of P_gs you plugged in values for the E field. These should be scaled as well. That brings the huge diagonal element for it down. After fixing this I get a condition number of 510. Which is perfectly reasonable.

 

[Ulver48]

In my code, I replaced all the Ereal, Eimag terms with 1e18*Ereal, 1e18*Ereal. After that, I computed the Jacobian Matrix ( which up to this point is wrong) and then I divided every element in the last two rows by 1e18. What have I done wrong ?

Pwl=0.271205495697359;
Pgs=0.500000000000001;
Pes=0.714813705790668;
Ereal=6.815768162566359e+17*1e-18   ;
Eimag=-4.788674222120425e+17*1e-18;

Delta=1.884955592153876e+10;
%%%%%%%%%%%%%%%%%%HERE COMES THE CODE%%%%%%%%%%%%%%%%%%%%%
% Fwl - Terms
dFwldPwl = -4.0833e9+Pes*1.8625e9;
dFwldPgs=0;
dFwldPes=2.2071e8+Pwl*1.8625e9;
dFwldEreal = 0;
dFwldEimag=0;

%Fgs - Terms
dFgsdPwl=0;
dFgsdPgs=-1.9514e-10*2*((Ereal*1e18)^2+(Eimag*1e18)^2)-3.75e11*Pes-8.3333e8-5.5556e10*(1-Pes);
dFgsdPes=3.75e11*(1-Pgs)+5.5556e10*Pgs;
dFgsdEreal=-1.9514e-10*(2*Ereal*1e18)*(2*Pgs-1);
dFgsdEimag=-1.9514e-10*(2*Eimag*1e18)*(2*Pgs-1);

%Fes - Terms
dFesdPwl=9.439e11*(1-Pes)+1e11*Pes;
dFesdPgs=1.6667e11*(1-Pes)+1.25e11*Pes;
dFesdPes=-9.439e11*Pwl-1.6667e11*Pgs-1e11*(1-Pwl)-1.25e11*(1-Pgs)-2e9;
dFesdEreal=0;
dFesdEimag=0;

%F_Ereal - Terms
dF_ErealdPwl=0;
dF_ErealdPgs=0.5*6.097e10*2*(Ereal*1e18-1.3*Eimag*1e18);
dF_ErealdPes=0;
dF_ErealdEreal=0.5*(6.0979e10*(2*Pgs-1)-6.3104e10);
dF_ErealdEimag=0.5*(6.0979e10*(2*Pgs-1)-6.3104e10)*(-1.3)-Delta;

%F_Eimag - Terms
dF_EimagdPwl=0;
dF_EimagdPgs=0.5*6.097e10*2*(Eimag*1e18+1.3*Ereal*1e18);
dF_EimagdPes=0;
dF_EimagdEreal=0.5*(6.0979e10*(2*Pgs-1)-6.3104e10)*(1.3)+Delta;
dF_EimagdEimag=0.5*(6.0979e10*(2*Pgs-1)-6.3104e10);

J=[dFwldPwl dFwldPgs dFwldPes dFwldEreal dFwldEimag;
    dFgsdPwl dFgsdPgs dFgsdPes dFgsdEreal dFgsdEimag;
    dFesdPwl dFesdPgs dFesdPes  dFesdEreal dFesdEimag;
    dF_ErealdPwl/1e18 dF_ErealdPgs/1e18 dF_ErealdPes/1e18 dF_ErealdEreal/1e18 dF_ErealdEimag/1e18;
    dF_EimagdPwl/1e18 dF_EimagdPgs/1e18 dF_EimagdPes/1e18 dF_EimagdEreal/1e18 dF_EimagdEimag/1e18];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Edit: I can see some mistakes there. I forgot adding the 1e18 term in the dF_EimagdEreal, dF_EimagdEimag tems and so on. I will fix it and then post the result. My bad.

 

[Ulver48]

Ok, I fixed my mistakes. Now we get the same matrix except from the second row. How did you compute it ?

Edit: I see. There is a second mistake. Oh my god. If only I wasn't so hasty with my math.

 

[Ulver48]

Ok. My only problem is the element in the second row and second column. You get -2.7080e-10, but shouldn't it be -2.7080e26 ?

 

[mfb]

The previous value was 10²⁶ and we scaled it by the square of a factor 10¹⁸. 26-18-18 = -10.

 

[Ulver48]

Yes but we now use the values $E^n_{real}=0.6815$, $E^n_{imag}=-0.4789$. So this term now is
$$ -1.9514*10^{-10}*2*((E^n_{real}*10^{18})^2+(E^n_{imag}*10^{18})^2)+\ldots$$
The exponent 26 is there again. 18+18-10=36-10=26.

 

[mfb]

Okay, if that is true your system is probably poorly conditioned. P_gs will go to 1/2 on a timescale much shorter than all other changes. For all practical purposes it will stay at 1/2 while the others change afterwards. It might be possible to plug that value in as a constant.

 

[Ulver48]

Yes, I think you are right. Solving the set of ODEs with a stiff solver in Matlab, it seems that the $P_{gs}$ reaches the value 0.5 sooner than the other variables reach their steady state. The problem is that I cannot ommit the ODE associated with the $P_{gs}$ variable. In order to create my bifurcation diagram I will be changing the parameters $\Delta\omega, k_{eff}$ and they will affect the variables $E_{rea},E_{imag}$ , which are terms inside the $F_{gs}$ function.