Scipy.optimize.root does not converge in Python while Matlab fsolve works, why?

2024/9/22 15:24:52

I am trying to find the root y of a function called f using Python.

Here is my code:

def f(y):w,p1,p2,p3,p4,p5,p6 = y[:7] t1 = w - 0.99006633*(p1**0.5) - (-1.010067)*((1-p1))t2 = w - 22.7235687*(p2**0.5) - (-1.010067)*((1-p2))t3 = w - 9.71323491*(p3**0.5) - (-1.010067)*((1-p3))t4 = w - 2.43852877*(p4**0.5) - (-1.010067)*((1-p4))t5 = w - 3.93640207*(p5**0.5) - (-1.010067)*((1-p5))t6 = w - 9.22688144*(p6**0.5) - (-1.010067)*((1-p6))t7 = p1 + p2 + p3 + p4 + p5 + p6 - 1return [t1,t2,t3,t4,t5,t6,t7]x0 = np.array([-0.01,0.4,0.1,0.2,0.1,0.1,0.1])
sol = scipy.optimize.root(f, x0)
print sol 

Python does not find the root. However there is one, I found it with the function fsolve in Matlab.

It is:

[ 0.3901, 0.6166, 0.0038, 0.0202, 0.2295, 0.1076, 0.0223]

I really want to use Python. Can anyone explain why scipy.optimize.root in Python does not converge while fsolve in Matlab does?

For info, scipy.optimize.solve does not converge either.

Answer

Try a different method. For me, method="lm" (I'm guessing Levenberg-Marquardt, but I'm not entirely sure) works very well:

import numpy as np
import scipy.optimizedef f(y):w,p1,p2,p3,p4,p5,p6 = y[:7]t1 = w - 0.99006633*(p1**0.5) - (-1.010067)*((1-p1))t2 = w - 22.7235687*(p2**0.5) - (-1.010067)*((1-p2))t3 = w - 9.71323491*(p3**0.5) - (-1.010067)*((1-p3))t4 = w - 2.43852877*(p4**0.5) - (-1.010067)*((1-p4))t5 = w - 3.93640207*(p5**0.5) - (-1.010067)*((1-p5))t6 = w - 9.22688144*(p6**0.5) - (-1.010067)*((1-p6))t7 = p1 + p2 + p3 + p4 + p5 + p6 - 1return [t1,t2,t3,t4,t5,t6,t7]x0 = np.array([-0.01,0.4,0.1,0.2,0.1,0.1,0.1])
sol = scipy.optimize.root(f, x0, method='lm')assert sol['success']
print 'Solution: ', sol.x
print 'Misfit: ', f(sol.x)

This yields:

Solution: [ 0.39012036  0.61656436  0.00377616  0.02017937  0.22954825 0.10763827  0.02229359]
Misfit: [0.0, 0.0, 1.1102230246251565e-16, -1.1102230246251565e-16,   1.1102230246251565e-16, 0.0, -2.2204460492503131e-16]

I'm actually a bit surprised Levenberg-Marquardt isn't the default. It's usually one of the first "gradient-descent" style solvers one would try.

https://en.xdnf.cn/q/71930.html

Related Q&A

Query CPU ID from Python?

How I can find processor id with py2.6, windows OS?I know that there is pycpuid, but I cant compile this under 2.6.

Recover from segfault in Python

I have a few functions in my code that are randomly causing SegmentationFault error. Ive identified them by enabling the faulthandler. Im a bit stuck and have no idea how to reliably eliminate this pro…

python and using self in methods

From what I read/understand, the self parameter is similiar to this.Is that true?If its optional, what would you do if self wasnt passed into the method?

How to extract only characters from image?

I have this type of image from that I only want to extract the characters.After binarization, I am getting this image img = cv2.imread(the_image.jpg) gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) thresh…

connect to Azure SQL database via pyodbc

I use pyodbc to connect to my local SQL database which works withoout problems.SQLSERVERLOCAL=Driver={SQL Server Native Client 11.0};Server=(localdb)\\v11.0;integrated security = true;DATABASE=eodba; c…

Generator function for prime numbers [duplicate]

This question already has answers here:Simple prime number generator in Python(28 answers)Closed 9 years ago.Im trying to write a generator function for printing prime numbers as followsdef getPrimes(n…

`ValueError: too many values to unpack (expected 4)` with `scipy.stats.linregress`

I know that this error message (ValueError: too many values to unpack (expected 4)) appears when more variables are set to values than a function returns. scipy.stats.linregress returns 5 values accord…

Django prefetch_related from foreignkey with manytomanyfield not working

For example, in Django 1.8:class A(models.Model):x = models.BooleanField(default=True)class B(models.Model):y = models.ManyToManyField(A)class C(models.Model):z = models.ForeignKey(A)In this scenario, …

Running SimpleXMLRPCServer in separate thread and shutting down

I have a class that I wish to test via SimpleXMLRPCServer in python. The way I have my unit test set up is that I create a new thread, and start SimpleXMLRPCServer in that. Then I run all the test, and…

Clean Python multiprocess termination dependant on an exit flag

I am attempting to create a program using multiple processes and I would like to cleanly terminate all the spawned processes if errors occur. below Ive wrote out some pseudo type code for what I think …