How to solve the problem without using symbolic computation

Question

I have the following simple nonlinear equations with two unknowns only:

$$\left\{ \begin{array}{c} \int_1^2{\dfrac{ e^{{a_1} x+{a_2} x^3}}{1+x^2}} \, dx=1 \\[13pt] \int_1^2{ x^2 e^{{a_1} x+{a_2} x^3}} \, dx=\dfrac{1}3 \\ \end{array} \right.$$

I understand it is almost impossible to obtain all the possible numerical solutions in the complex field $\mathbb{C}$. So any numerical complex solution per a specific given initial value with good enough accuracy would be OK.

By symbolic-numerical computation, it is easy to obtain one solution:

$$\left\{\quad \begin{array}{ccrcr} {a_1}&= & 2.083619981922075 &-& 2.163052112144277\; i \\ {a_2}&= & -0.144311264409679 &+& 1.206590594556891\; i \\ \end{array} \right.$$

However, I found it difficult to handle such a problem using the built-in integral functions in Matlab since such integral and quadgk functions as accept function handle-form only integrands.

My try is first composing the nonlinear equations into user defined functions, then solve it per iteration algorithm (e.g. Newton's method) with a specific initial value), but when I tried to pass symbolic variables or global variables $a_1$, $a_2$, problems (error messages) would happen.

Though symbolic computation is possible to obtain a solution, it is too slow especially for more general nonlinear equations. Since the solution I need is only numerical one, I wonder:

Is it possible to solve such a problem by Matlab without using the symbolic computation? If the integral function has to accept symbolics $a_1$, $a_2$, how to handle it in Matlab (or C++)?

Answer 1

You can solve this numerically in Python without symbolic computation.

from __future__ import print_function, division
import numpy as np
from numpy import exp
from scipy.integrate import quad
from scipy.optimize import root

def f1(a1, a2, x):
    return exp(a1 * x + a2 * x * x * x) / (1 + x * x)

def f2(a1, a2, x):
    return exp(a1 * x + a2 * x * x * x) * x

def g(a1, a2, f):
    def r(x):
        return f(a1, a2, x).real
    def i(x):
        return f(a1, a2, x).imag
    y = quad(r, 1, 2)[0] + quad(i, 1, 2)[0] * 1j
    return y

def func(Z):
    a1 = Z[0] + Z[1] * 1j
    a2 = Z[2] + Z[3] * 1j
    y1 = g(a1, a2, f1) - 1
    y2 = g(a1, a2, f2) - 1/3
    return np.array([y1.real, y1.imag, y2.real, y2.imag])

X0 = np.array([2.0, -2.0, -0.1, 1.0])
print(root(func, X0, tol=1e-12))

output:

fjac: array([[-0.30890916,  0.05007557,  0.09254457,  0.94525291],
       [-0.12460792, -0.41710616, -0.89760745,  0.06925448],
       [ 0.41309815, -0.83196711,  0.34050112,  0.14573825],
       [ 0.84758358,  0.36241635, -0.26418811,  0.28365668]])
     fun: array([ -2.88657986e-15,  -2.49800181e-15,  -1.89293026e-14,
        -2.43138842e-14])
 message: 'The solution converged.'
    nfev: 24
     qtf: array([ -9.32740208e-12,   6.42466839e-12,  -3.53520661e-12,
        -1.78427395e-12])
       r: array([ -3.43272029,   0.72043771, -14.00467579, -11.98268077,
        -2.65005807,  13.29034808, -12.53526771,  -7.04373163,
         2.58188876,  -7.03891799])
  status: 1
 success: True
       x: array([ 2.08361998, -2.16305211, -0.14431126,  1.20659059])

Notice that I've changed the definition of the second function to use an x term instead of an x^2 term so that it matches your solution :) I guess this was a typo in your question.