c语言栈求解表达式_非线性方程组的编程求解方法

【作者声明】

　　本文所有文字均为作者原创，所有图片均为作者本人亲自拍摄或制作。

前言

　　我们在以前曾经写过一篇文章，是关于科学计算器fx-991CN X求解非线性方程组的内容，主要的原理是把一个二元非线性方程组的每一个表达式（可以看作是隐函数）进行显化，将某个未知数变量使用一个函数来表示，代入另一个方程，求解出一个未知数之后，再代回原来的方程，求解出另一个未知数。

電卓院亜紀良：fx-991CN X非线性方程组的求解应用zhuanlan.zhihu.com

　　本次我们以一个形式上比较简单的非线性方程组

$\left\{\begin{matrix} x^3+y^3-4=0\\ x^4+y^2-3=0 \end{matrix}\right.$ 为例进行讲解。

　　能够直接求解二元非线性方程组的计算器当然是有的，不过都是带有CAS功能的图形计算器。例如TI-Nspire CX CAS，其SOLVE指令可以直接计算，也可以通过图形的方法来求解。

　　不过，我们这里要讲的当然不是这种方法，毕竟带有CAS功能的计算器的用户不多，大多数图形计算器或可编程计算器都没有CAS功能，然而二元非线性方程组又是经常需要求解的方程类型之一，所以我们考虑通过编程的方式，来实现对任意的二元非线性方程组的求解。

从图形解法到程序解法

　　一般的图形计算器如果要采用图形解法，也必须把方程进行显化。方程的显化难度视方程的复杂程度有所不同，如果遇到显化比较麻烦的方程，这种方法也比较吃力。例如示例中的方程，将其显化之后，就要在图形计算器上输入三个表达式：

$\left\{\begin{matrix} y_1=\sqrt[3]{4-x^3}\\ y_2=\sqrt{3-x^4}\\ y_3=-\sqrt{3-x^{4}} \end{matrix}\right.$

然后对这三个表达式进行绘图，再求交点，就像这样：

CASIO fx-CG50

TI-83/84 Plus

　　图形方法虽然直观，但仍然需要对方程进行显化，如果能够按照方程组的形式直接输入，求解的效率就大幅提高了。我们知道，fx-991CN X的SOLVE功能是使用牛顿法求解的。实际上对于二元非线性方程组，同样可以使用牛顿法进行求解。

【求解原理】二元非线性方程组的牛顿迭代法

　　设有非线性方程组

$\left\{\begin{matrix} u(x,y)=0\\ v(x,y)=0 \end{matrix}\right.$ ，记雅可比行列式

$J=\begin{vmatrix} \frac{\partial u}{\partial x} & \frac{\partial u}{\partial y}\\ \frac{\partial v}{\partial x} & \frac{\partial v}{\partial y} \end{vmatrix}$ ，那么有

$x_{n+1}=x_n+\frac{1}{J}\begin{vmatrix} \frac{\partial u}{\partial y} & u_n\\ \frac{\partial v}{\partial y} & v_n \end{vmatrix}$ ，

$y_{n+1}=y_n+\frac{1}{J}\begin{vmatrix} u_n & \frac{\partial u}{\partial x}\\ v_n & \frac{\partial v}{\partial x} \end{vmatrix}$ 。

　　在程序中实现的时候，偏导数用数值方法计算，这样就可以避免求导函数的计算。

　　例如计算

$\frac{\partial u}{\partial x}|_{x=x_0}$ ，令

$y=y_0$ ，然后计算

$\frac{u(x_0+h,y_0)-u(x_0,y_0)}{h}$ 即可，其中

$h$ 取一个很小的正数。

程序求解算法

　　为保护求解程序不被修改，仅作为调用入口使用，我们将求解程序和方程组分别用两个程序文件编写，其中求解程序命名为“SOLVE2”，方程组程序命名为“EQN2”。不同的计算器程序语法略有不同，此处提供通用的版本（接近CASIO计算器的类BASIC语言），实际编写时需要根据使用的计算器修改。

求解程序【SOLVE2】
"X0="?→A:"Y0="?→B:1E-5→S（指定求解的初始值，以及求导时的增量）
If A=0:Then S→A:IfEnd
If B=0:Then S→B:IfEnd（为解决第一次迭代时指定初始值(0, 0)让迭代条件中的A和B报错）
Lbl 0（设置标签）
A→X:B→Y
Prog "EQN2"（计算u和v）
P→G:Q→H
A+S→X:B→Y
Prog "EQN2"
(P-G)÷S→I:(Q-H)÷S→J（计算u和v对x的偏导数）
A→X:B+S→Y
Prog "EQN2"
(P-G)÷S→K:(Q-H)÷S→L（计算u和v对y的偏导数）
If Abs((IL-JK)^(-1)×(KH-GL)÷A)<1E-10 And Abs((IL-JK)^(-1)×(GJ-HI)÷B)<1E-10（迭代终止条件）
Then
"X=":A◢
"Y=":B◢
"END"
Else
A+(IL-JK)^(-1)×(KH-GL)→A（迭代u）
B+(IL-JK)^(-1)×(GJ-HI)→B（迭代v）
Goto 0（返回标签）
IfEnd方程组程序【EQN2】
X^3+Y^3-4→P（计算函数u(x,y)的值）
X^4+Y^2-3→Q（计算函数v(x,y)的值）

程序运行效果

　　这一程序每次使用时仅需要改变EQN2里面的程序内容，即修改方程，但要注意方程的形式一定要是等号的一侧是0，然后只需要输入方程另一侧的内容。运行程序时，为x和y赋值初始迭代值，等待求解即可。如果方程组有多个解，可以尝试多次运行程序，为x和y赋不同的初始值，得到不同的解。

　　例如指定初始值(1, 1)，得到的解为x=0.9676，y=1.4572，与图形解法一致。

CASIO fx-CG50

TI-83/84 Plus

CASIO fx-5800P

总结

　　二元非线性方程组的求解是经常遇到的一类问题之一，使用可编程计算器或者图形计算器，运用牛顿法的原理，即可编写出实用的求解程序，不再需要进行方程的显化计算，提高了解决理论或实际问题的效率。这一实例再一次体现可编程计算器的潜力之巨大。