## Coupled nonlinear equations

| categories: nonlinear algebra | tags: | View Comments

Suppose we seek the solution to this set of equations:

\begin{align} y &=& x^2 \\ y &=& 8 - x^2 \end{align}

To solve this we need to setup a function that is equal to zero at the solution. We have two equations, so our function must return two values. There are two variables, so the argument to our function will be an array of values.

from scipy.optimize import fsolve

def objective(X):
x, y = X            # unpack the array in the argument
z1 = y - x**2       # first equation
z2 = y - 8 + x**2   # second equation
return [z1, z2]     # list of zeros

x0, y0 = 1, 1           # initial guesses
guess = [x0, y0]
sol = fsolve(objective, guess)
print sol

# of course there may be more than one solution
x0, y0 = -1, -1           # initial guesses
guess = [x0, y0]
sol = fsolve(objective, guess)
print sol

[ 2.  4.]
[-2.  4.]


org-mode source

## Finding the nth root of a periodic function

| categories: nonlinear algebra | tags: heat transfer | View Comments

There is a heat transfer problem where one needs to find the n^th root of the following equation: $$x J_1(x) - Bi J_0(x)=0$$ where $$J_0$$ and $$J_1$$ are the Bessel functions of zero and first order, and $$Bi$$ is the Biot number. We examine an approach to finding these roots.

First, we plot the function.

from scipy.special import jn, jn_zeros
import matplotlib.pyplot as plt
import numpy as np

Bi = 1

def f(x):
return x * jn(1, x) - Bi * jn(0, x)

X = np.linspace(0, 30, 200)
plt.plot(X, f(X))
plt.savefig('images/heat-transfer-roots-1.png')


You can see there are many roots to this equation, and we want to be sure we get the n^{th} root. This function is pretty well behaved, so if you make a good guess about the solution you will get an answer, but if you make a bad guess, you may get the wrong root. We examine next a way to do it without guessing the solution. What we want is the solution to $$f(x) = 0$$, but we want all the solutions in a given interval. We derive a new equation, $$f'(x) = 0$$, with initial condition $$f(0) = f0$$, and integrate the ODE with an event function that identifies all zeros of $$f$$ for us. The derivative of our function is $$df/dx = d/dx(x J_1(x)) - Bi J'_0(x)$$. It is known (http://www.markrobrien.com/besselfunct.pdf) that $$d/dx(x J_1(x)) = x J_0(x)$$, and $$J'_0(x) = -J_1(x)$$. All we have to do now is set up the problem and run it.

from pycse import *  # contains the ode integrator with events

from scipy.special import jn, jn_zeros
import matplotlib.pyplot as plt
import numpy as np

Bi = 1

def f(x):
"function we want roots for"
return x * jn(1, x) - Bi * jn(0, x)

def fprime(f, x):
"df/dx"
return x * jn(0, x) - Bi * (-jn(1, x))

def e1(f, x):
"event function to find zeros of f"
isterminal = False
value = f
direction = 0
return value, isterminal, direction

f0 = f(0)
xspan = np.linspace(0, 30, 200)

x, fsol, XE, FE, IE = odelay(fprime, f0, xspan, events=[e1])

plt.plot(x, fsol, '.-', label='Numerical solution')
plt.plot(xspan, f(xspan), '--', label='Analytical function')
plt.plot(XE, FE, 'ro', label='roots')
plt.legend(loc='best')
plt.savefig('images/heat-transfer-roots-2.png')

for i, root in enumerate(XE):
print 'root {0} is at {1}'.format(i, root)

plt.show()

root 0 is at 1.25578377377
root 1 is at 4.07947743741
root 2 is at 7.15579904465
root 3 is at 10.2709851256
root 4 is at 13.3983973869
root 5 is at 16.5311587137
root 6 is at 19.6667276775
root 7 is at 22.8039503455
root 8 is at 25.9422288192
root 9 is at 29.081221492


You can work this out once, and then you have all the roots in the interval and you can select the one you want.

org-mode source

## Method of continuity for solving nonlinear equations - Part II

| categories: nonlinear algebra | tags: | View Comments

Matlab post Yesterday in Post 1324 we looked at a way to solve nonlinear equations that takes away some of the burden of initial guess generation. The idea was to reformulate the equations with a new variable $$\lambda$$, so that at $$\lambda=0$$ we have a simpler problem we know how to solve, and at $$\lambda=1$$ we have the original set of equations. Then, we derive a set of ODEs on how the solution changes with $$\lambda$$, and solve them.

Today we look at a simpler example and explain a little more about what is going on. Consider the equation: $$f(x) = x^2 - 5x + 6 = 0$$, which has two roots, $$x=2$$ and $$x=3$$. We will use the method of continuity to solve this equation to illustrate a few ideas. First, we introduce a new variable $$\lambda$$ as: $$f(x; \lambda) = 0$$. For example, we could write $$f(x;\lambda) = \lambda x^2 - 5x + 6 = 0$$. Now, when $$\lambda=0$$, we hve the simpler equation $$- 5x + 6 = 0$$, with the solution $$x=6/5$$. The question now is, how does $$x$$ change as $$\lambda$$ changes? We get that from the total derivative of how $$f(x,\lambda)$$ changes with $$\lambda$$. The total derivative is:

$$\frac{df}{d\lambda} = \frac{\partial f}{\partial \lambda} + \frac{\partial f}{\partial x}\frac{\partial x}{\partial \lambda}=0$$

We can calculate two of those quantities: $$\frac{\partial f}{\partial \lambda}$$ and $$\frac{\partial f}{\partial x}$$ analytically from our equation and solve for $$\frac{\partial x}{\partial \lambda}$$ as

$$\frac{\partial x}{\partial \lambda} = -\frac{\partial f}{\partial \lambda}/\frac{\partial f}{\partial x}$$

That defines an ordinary differential equation that we can solve by integrating from $$\lambda=0$$ where we know the solution to $$\lambda=1$$ which is the solution to the real problem. For this problem: $$\frac{\partial f}{\partial \lambda}=x^2$$ and $$\frac{\partial f}{\partial x}=-5 + 2\lambda x$$.

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def dxdL(x, Lambda):
return -x**2 / (-5.0 + 2 * Lambda * x)

x0 = 6.0/5.0
Lspan = np.linspace(0, 1)
x = odeint(dxdL, x0, Lspan)

plt.plot(Lspan, x)
plt.xlabel('$\lambda$')
plt.ylabel('x')
plt.savefig('images/nonlin-contin-II-1.png')


We found one solution at x=2. What about the other solution? To get that we have to introduce $$\lambda$$ into the equations in another way. We could try: $$f(x;\lambda) = x^2 + \lambda(-5x + 6)$$, but this leads to an ODE that is singular at the initial starting point. Another approach is $$f(x;\lambda) = x^2 + 6 + \lambda(-5x)$$, but now the solution at $$\lambda=0$$ is imaginary, and we do not have a way to integrate that! What we can do instead is add and subtract a number like this: $$f(x;\lambda) = x^2 - 4 + \lambda(-5x + 6 + 4)$$. Now at $$\lambda=0$$, we have a simple equation with roots at $$\pm 2$$, and we already know that $$x=2$$ is a solution. So, we create our ODE on $$dx/d\lambda$$ with initial condition $$x(0) = -2$$.

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def dxdL(x, Lambda):
return (5 * x - 10) / (2 * x - 5 * Lambda)

x0 = -2
Lspan = np.linspace(0, 1)
x = odeint(dxdL, x0, Lspan)

plt.plot(Lspan, x)
plt.xlabel('$\lambda$')
plt.ylabel('x')
plt.savefig('images/nonlin-contin-II-2.png')


Now we have the other solution. Note if you choose the other root, $$x=2$$, you find that 2 is a root, and learn nothing new. You could choose other values to add, e.g., if you chose to add and subtract 16, then you would find that one starting point leads to one root, and the other starting point leads to the other root. This method does not solve all problems associated with nonlinear root solving, namely, how many roots are there, and which one is “best” or physically reasonable? But it does give a way to solve an equation where you have no idea what an initial guess should be. You can see, however, that just like you can get different answers from different initial guesses, here you can get different answers by setting up the equations differently.

org-mode source

## Finding equilibrium conversion

| categories: nonlinear algebra | tags: | View Comments

A common problem to solve in reaction engineering is finding the equilibrium conversion.1 A typical problem to solve is the following nonlinear equation:

$$1.44 = \frac{X_e^2}{(1-X_e)^2}$$

To solve this we create a function:

$$f(X_e)=0=1.44 - \frac{X_e^2}{(1-X_e)^2}$$

and use a nonlinear solver to find the value of $$X_e$$ that makes this function equal to zero. We have to provide an initial guess. Chemical intuition suggests that the solution must be between 0 and 1, and mathematical intuition suggests the solution might be near 0.5 (which would give a ratio near 1).

Here is our solution.

from scipy.optimize import fsolve

def func(Xe):
z = 1.44 - (Xe**2)/(1-Xe)**2
return z

X0 = 0.5
Xe, = fsolve(func, X0)
print('The equilibrium conversion is X = {0:1.2f}'.format(Xe))

The equilibrium conversion is X = 0.55


org-mode source

## Counting roots

| categories: nonlinear algebra | tags: | View Comments

Matlab post The goal here is to determine how many roots there are in a nonlinear function we are interested in solving. For this example, we use a cubic polynomial because we know there are three roots.

$$f(x) = x^3 + 6x^2 - 4x -24$$

## 1 Use roots for this polynomial

This ony works for a polynomial, it does not work for any other nonlinear function.

import numpy as np
print np.roots([1, 6, -4, -24])

[-6.  2. -2.]


Let us plot the function to see where the roots are.

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(-8, 4)
y = x**3 + 6 * x**2 - 4*x - 24
plt.plot(x, y)
plt.savefig('images/count-roots-1.png')


Now we consider several approaches to counting the number of roots in this interval. Visually it is pretty easy, you just look for where the function crosses zero. Computationally, it is tricker.

## 2 method 1

Count the number of times the sign changes in the interval. What we have to do is multiply neighboring elements together, and look for negative values. That indicates a sign change. For example the product of two positive or negative numbers is a positive number. You only get a negative number from the product of a positive and negative number, which means the sign changed.

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(-8, 4)
y = x**3 + 6 * x**2 - 4*x - 24

print np.sum(y[0:-2] * y[1:-1] < 0)

3


This method gives us the number of roots, but not where the roots are.

## 3 Method 2

Using events in an ODE solver python can identify events in the solution to an ODE, for example, when a function has a certain value, e.g. f(x) = 0. We can take advantage of this to find the roots and number of roots in this case. We take the derivative of our function, and integrate it from an initial starting point, and define an event function that counts zeros.

$$f'(x) = 3x^2 + 12x - 4$$

with f(-8) = -120

import numpy as np
from pycse import odelay

def fprime(f, x):
return 3.0 * x**2 + 12.0*x - 4.0

def event(f, x):
value = f # we want f = 0
isterminal = False
direction = 0
return value, isterminal, direction

xspan = np.linspace(-8, 4)
f0 = -120

X, F, TE, YE, IE = odelay(fprime, f0, xspan, events=[event])
for te, ye in zip(TE, YE):
print 'root found at x = {0: 1.3f}, f={1: 1.3f}'.format(te, ye)

root found at x = -6.000, f=-0.000
root found at x = -2.000, f=-0.000
root found at x =  2.000, f= 0.000