-1

I'm trying to solve a system of 30 non linear equations (unknown variables sit inside exponential functions). I have exactly 30 variables and 30 unknowns. I'm trying to use scipy.optimize.fsolve but it is essentially not doing anything. I'll post my attempt below.

Below, yu,yd,and Vckm are known matrices with elements ranging from 1e-6 to 1, and the function qyuk(ca,cb,lam)depends on Exp[(1 - ca - cb)lam]. My system is complicated because this exponential has to reproduce small numbers (~1e-6) so my solution is extremely sensitive. 

def f(c):
    cq1,cq2,cq3,cu1,cu2,cu3,cd1,cd2,cd3,ce1,ce2,ce3,lu1,lu2,lu3,lu4,lu5,lu6,lu7,lu8,lu9,du1,du2,du3,du4,du5,du6,du7,du8,du9 = c


    Gu = [[qyuk(cq1,cu1,lam),qyuk(cq1,cu2,lam),qyuk(cq1,cu3,lam)],
          [qyuk(cq2,cu1,lam),qyuk(cq2,cu2,lam),qyuk(cq2,cu3,lam)],
          [qyuk(cq3,cu1,lam),qyuk(cq3,cu2,lam),qyuk(cq3,cu3,lam)]]

    Gd = [[qyuk(cq1,cd1,lam),qyuk(cq1,cd2,lam),qyuk(cq1,cd3,lam)],
          [qyuk(cq2,cd1,lam),qyuk(cq2,cd2,lam),qyuk(cq2,cd3,lam)],
          [qyuk(cq3,cd1,lam),qyuk(cq3,cd2,lam),qyuk(cq3,cd3,lam)]]

    Gu_squared = np.matmul(Gu,np.transpose(Gu))
    Gd_squared = np.matmul(Gd,np.transpose(Gd))

    Ul = [[lu1,lu2,lu3],[lu4,lu5,lu6],[lu7,lu8,lu9]]
    Dl = [[du1,du2,du3],[du4,du5,du6],[du7,du8,du9]]

    rul = np.matmul(np.transpose(Ul),np.matmul(yu,np.matmul(yu,Ul)))
    rdl = np.matmul(np.transpose(Dl),np.matmul(yd,np.matmul(yd,Dl)))

    rvckm = np.matmul(np.transpose(Ul),Dl)


    return (
        Gu_squared[0][0] - rul[0][0],Gu_squared[0][1] - rul[0][1],Gu_squared[0][2] - rul[0][2],
        Gu_squared[1][0] - rul[1][0],Gu_squared[1][1] - rul[1][1],Gu_squared[1][2] - rul[1][2],
        Gu_squared[2][0] - rul[2][0],Gu_squared[2][1] - rul[2][1],Gu_squared[2][2] - rul[2][2],
        Gd_squared[0][0] - rdl[0][0],Gd_squared[0][1] - rdl[0][1],Gd_squared[0][2] - rdl[0][2],
        Gd_squared[1][0] - rdl[1][0],Gd_squared[1][1] - rdl[1][1],Gd_squared[1][2] - rdl[1][2],
        Gd_squared[2][0] - rdl[2][0],Gd_squared[2][1] - rdl[2][1],Gd_squared[2][2] - rdl[2][2],
        qyuk(ce1,ce1,lam) - ye[0][0],qyuk(ce2,ce2,lam) - ye[1][1],qyuk(ce3,ce3,lam) - ye[2][2],
        Vckm[0][0] - rvckm[0][0],Vckm[0][1] - rvckm[0][1],Vckm[0][2] - rvckm[0][2],
        Vckm[1][0] - rvckm[1][0],Vckm[1][1] - rvckm[1][1],Vckm[1][2] - rvckm[1][2],
        Vckm[2][0] - rvckm[2][0],Vckm[2][1] - rvckm[2][1],Vckm[2][2] - rvckm[2][2],
        )

c = scipy.optimize.fsolve(f,np.full((30,1),.6))

My problem is that fsolve is returning garbage and my solution doesn't solve my system of equations, not even close. Any help would be appreciated. Thank you.

petezurich
  • 9,280
  • 9
  • 43
  • 57
CStarAlgebra
  • 197
  • 2
  • 10
  • 1
    30 equations is not large. It's routine to solve problem with millions of degrees of freedom. Non-linear problems may not have a solution; you may choose a bad initial starting position; your matricies may be ill-conditioned. Can't answer based on what you've posted. – duffymo Sep 28 '18 at 17:50

1 Answers1

0

Since your system is very sensitive to its parameters and coefficients, the starting point for an iterative process is very important. If that point was selected arbitrarily, the Newton-Raphson iterations may not converge at all. Besides, non-linear systems may have no solutions at all, of have more than one solution. If you can somehow provide the full system with coefficients, it will be possible to begin working with that.

SparseSolverCodes
  • 61
  • 3