0

I am a researcher in the field of solid mechanics and new to this website as far as asking the question side is concerned.

I am trying to use nsolve of sympy to find solutions of system of non-linear algebraic equations. My system has three non-linear equations with three unknowns. Each of these equations also has an input parameter 'V' which I specify explicitly.

i = 7.0
NLE_1 = Nonlinear_eq_1.subs(V, i)
NLE_2 = Nonlinear_eq_2.subs(V, i)
NLE_3 = Nonlinear_eq_3.subs(V, i)

As said earlier, Nonlinear_eq_1, Nonlinear_eq_2 and Nonlinear_eq_3 are functions of a1, a2 and a3. Now this non-linear system has upper threshold value of V. I need to find out this critical V value and which is the actual unknown that I am interested in.

For given system parameters, I have tried to zoom in on the critical V value as follows:

In the beginning, I iterated i manually (without employing any looping e.g., for loop) so that it increases from 0 with increment of 1. Using the following line of code, I found out the roots of the system of non-linear equations:

ans_nonlinear = smp.nsolve([NLE_1, NLE_2, NLE_3], [a1, a2, a3], linear_answers)

where, 'linear_answers' is a tuple which is obtained by solving the Nonlinear_eq_1, Nonlinear_eq_2 and Nonlinear_eq_3 by only considering linear terms involving a1, a2, a3 as follows:

ans_linear = smp.solve((LE_1, LE_2, LE_3), (a1, a2, a3))
linear_answers = (float(ans_linear.get(a1)), float(ans_linear.get(a2)), float(ans_linear.get(a3)))

Now, from V = 0 till V = 7, I got definite answers for ans_nonlinear. However, when V = 8, nsolve gave the following error:

Could not find root within given tolerance. (2.16855306310932770901e-19 > 2.16840434497100886801e-19) Try another starting point or tweak arguments.

This error is indirectly stating that V has a critical value inbetween 7 and 8.

I wanted to know how to use this error so that I can raise a user-defined flag and suppress the error as well. I am intending on using this flag so that I can then go back to V = 7 and start incrementing V as 7.1, 7.2, 7.3 etc. Using this logic in a loop, I will then hone in on the critical value of V to the desired accuracy.

This problem is basically on finding out static pull-in instability parameters of micro-cantilevers using three-mode reduced order model (ROM).

James Z
  • 12,209
  • 10
  • 24
  • 44
iKD
  • 1
  • 1
  • You could catch the exception using `try/except`. Your approach does not sound very robust. I attempt the general problem differently e.g. by searching for the minimum of the equations or something. Are your nonlinear equations polynomial or do they have transcendental functions? – Oscar Benjamin Mar 02 '23 at 20:54
  • Thank you, Oscar Benjamin, for your suggestion. Each of the three non-linear equations has 71 terms. These equations can not be classified as transcendental. However, each equation involves terms such as a1^4, a1^2*a2*a3, a1*a2*a3*a4 etc. I have tried searching for a minimum approach to find out critical V, but it is not appropriate as system is highly non-linear and involves multiple extremum values. Will you kindly let me know the way to use try/except in more detail with regard to my query? – iKD Mar 02 '23 at 21:05
  • This problem is not reproducible without your actual equations – Reinderien Mar 02 '23 at 23:19
  • Try/except is a standard Python feature. You haven't shown any code that can be adapted to add try/except but any beginner Python tutorial will show you how to use try/except. If your equations are polynomial in `a1`, `a2`, `a3` and also `V` then you can use techniques based on resultants and discriminants to identify the critical values. I would probably attempt that in the first instance. If it isn't possible to show your full system here it might be helpful to instead show complete code for a simplified version. – Oscar Benjamin Mar 03 '23 at 11:58
  • @Oscar Benjamin relevant section of my code is as follows: – iKD Mar 03 '23 at 12:49

1 Answers1

0

I would recommend that you use bisection to find the critical value. Your function for bisection could return a solution if found or None if it is not. Stop bisecting when your bounds are close enough together. e.g. if you had 1 equations eq that needed a parameter v specified and you had a guess g for the solution. Here is a toy 1-D example:

from sympy import nsolve
from sympy.abc import x, v

def f(vi, g):
    try: return nsolve(eq.subs(v, vi), g)
    except ValueError: return None

eq = x**2 + 1/10**v
vL = 100
vR = 0
g = 1  # some guess
while vR - vL > 1e-3:
    vm = (vL + vR)/2
    newg = f(vm, g)
    if newg is None:
        vR = vm
    else:
        g = newg
        vL = vm

print('v between',vL,'and',vR)

outputs

v between 9.803009033203125 and 9.80224609375
smichr
  • 16,948
  • 2
  • 27
  • 34