软件OpenLu中的解方程函数

 

  >>> [OpenLu的使用方法]

 

  欢迎加解方程\优化\拟合QQ群：562025053

 

    方程（组）的求解，难易程度差别较大。在OpenLu中，普通的方程（组）可借助LuMath库中的拟牛顿法netn和对分法btFindRoot求解，难度大的方程（组）须借助优化库LuOpt中的iFind、Find和Opt函数求解。

    （1）math::btFindRoot(f,a,b,h,k,eps)：对分法求非线性方程的实根

    f：自定义一元函数句柄。函数自变量不能重新赋值。
a,b：求根区间的左端点和右端点（b>a）。
h：搜索求根时采用的步长（h>0）。
k: 可能的实根个数。可缺省，缺省值为10。
eps：控制精度（eps>0）。可缺省，缺省值为1e-6。

    返回值：返回一个数组（存放搜索到的实根），或者返回nil。若返回数组的长度为k，则有可能在求根区间[a,b]内的实根未搜索完。

    例子：

f(x)=x^6-5*x^5+3*x^4+x^3-7*x^2+7*x-20;  //函数定义
math::btFindRoot[@f,-2.,5.,0.2].o[];   //函数o用于输出一个对象

    结果：

-1.40246 4.33376

    （2）math::netn(f,x,eps,t,h,k)：求非线性方程组一组实根的拟牛顿法

    f：自定义二元或2n(n>1)元函数句柄，由用户自编。

    如果f是二元函数，则两个参数为等长的一维数组：

f(x,y)=   //二元函数，x[i]为自变量，y[i]为方程右端点函数值
{
y[0]=f1(x[0],x[1],...,x[n-1]),
y[1]=f2(x[0],x[1],...,x[n-1]),
... ...
y[n-1]=fn(x[0],x[1],...,x[n-1])
};

    或者：

f(x1,x2,...,xn,y1,y2,...,yn)=   //2n元函数，xi为自变量，yi为方程右端点函数值
{
y1=f1(x1,x2,...,xn),
y2=f2(x1,x2,...,xn),
... ...
yn=fn(x1,x2,...,xn)
};

    x：双精度实型一维数组，长度不小于n，存放一组初值，返回方程组的一组实根。
eps：控制精度（eps>0）。可缺省该参数，缺省值eps=1e-6。
t：控制h大小的变量，0<t<1。可缺省该参数，缺省值h=0.1。
h：增量初值，在本函数中将被破坏。可缺省该参数，缺省值h=0.1。
k：允许的最大迭代次数。可缺省该参数，缺省值k=100。

    返回值：返回值为实际迭代次数。若返回值为-1表示代数方程奇异，若返回值为-2表示β=0，这两种情况可放宽精度要求、改变各个初值或改变各个方程顺序再试；若返回值等于0说明迭代了k次仍未满足精度要求，程序工作失败。

    例子：解下例方程组

x1*x1+x2*x2+x3*x3-1.0=0
2.0*x1*x1+x2*x2-4.0*x3=0
3.0*x1*x1-4.0*x2+x3*x3=0

    代码：

f(x1,x2,x3,y1,y2,y3)=  //函数定义
{
y1=x1*x1+x2*x2+x3*x3-1.0,
y2=2.0*x1*x1+x2*x2-4.0*x3,
y3=3.0*x1*x1-4.0*x2+x3*x3
};
!!!using["math"];      //使用命名空间math
main(:x,i)={
//申请一维数组并赋初值的常规方法。函数new表示申请一个对象，real_s表示数组对象，data后的数据为初值1.0,1.0,1.0
    x=new{real_s,data : 1.0,1.0,1.0},
    i=netn[@f,x],      //拟牛顿法解方程
    x.outa(),          //用函数outa输出一个数组
    i                  //返回迭代次数
};

    或者：

f(x,y)=
{
y[0]=x[0]*x[0]+x[1]*x[1]+x[2]*x[2]-1.0,
y[1]=2.0*x[0]*x[0]+x[1]*x[1]-4.0*x[2],
y[2]=3.0*x[0]*x[0]-4.0*x[1]+x[2]*x[2]
};
!!!using["math"];
main(:x,i)={
//申请一维数组并赋初值的简便方法。函数ra1(realarray1的缩写)表示申请一维数组，初值为1,1.0,1（ra1自动将整数转换为实数）
    i=netn[@f,x=ra1(1,1.0,1)],
x.o(),   //用函数o输出一个对象
    i
};

    结果：

0.785197 0.496611 0.369923

4

    （3）luopt::iFind(f, luopt::optmax,m, luopt::optrange,min,max, luopt::opteps,eps, luopt::optpara,x1,x2,..., luopt::optthis,i)：求单变量方程的全部解

    f：自定义n元函数句柄，不可缺省。格式如下：

f(x1,x2,...,xn)=
{
... ...
};

    默认求方程f(x1,x2,...,xn)=0第一个自变量的全部解，而其他自变量赋值为0.0。可以用参数optthis指定求解的自变量，也可以用参数optpara给出其他自变量的值。

    luopt::optmax,m：区间分割数目（大于等于10），区间分割数目越多精度越高。可以缺省该参数，缺省值为200。
luopt::optrange,min,max：指定求解区间。若缺省该参数，则min为-1e50，max为1e50。
luopt::opteps,eps：控制精度要求，eps>0。可以缺省该参数，缺省值为1e-6。
luopt::optpara,x1,x2,...：给指定求解自变量之外的其他自变量赋值，参数x1,x2,...的个数比全部自变量个数少1。若缺省该参数，其他自变量缺省值均为0.0。
luopt::optthis,i：指定求解的自变量。0表示第一个自变量；1表示第二个自变量，以此类推。若缺省该参数，i为0。

    返回值：解的个数。

    例子1：求方程的全部解：f(x)=x^6-5*x^5+3*x^4+x^3-7*x^2+7*x-20;

f(x)=x^6-5*x^5+3*x^4+x^3-7*x^2+7*x-20;
luopt::iFind[@f];

    结果（最后一个值为误差，下同）：

-1.402463030422577 3.552713678800501e-015

4.333755446919994 -2.984279490192421e-013

2

    例子2：求方程的全部解（x取-8～8，y取1）：f(x,y)=(x^2+y^2)^3-36*(x^2-y^2)^2;

!!!using("luopt");
f(x,y)=(x^2+y^2)^3-36*(x^2-y^2)^2;
iFind[@f,optrange,-8.,8.,optpara,1.];

    结果：

-5.530393535267971  1. -3.637978807091713e-012

-1.32936186218296   1. -7.105427357601002e-015

-0.8047014256478491 1. -1.776356839400251e-015

0.8047014256478491  1. -1.776356839400251e-015

1.32936186218296    1. -7.105427357601002e-015

5.530393535267971   1. -3.637978807091713e-012

6

    例子3：求方程的全部解（x取1，y取-8～8）：f(x,y)=(x^2+y^2)^3-36*(x^2-y^2)^2;

!!!using("luopt");
f(x,y)=(x^2+y^2)^3-36*(x^2-y^2)^2;
iFind[@f, optrange,-8.,8., optthis,1, optpara,1.];

    结果：

1. -5.530393535267971  -3.637978807091713e-012

1. -1.32936186218296   -7.105427357601002e-015

1. -0.8047014256478491 -1.776356839400251e-015

1. 0.8047014256478491  -1.776356839400251e-015

1. 1.32936186218296    -7.105427357601002e-015

1. 5.530393535267971   -3.637978807091713e-012

6

    例子4：怎么解出这个方程的所有解（很少有软件能顺利求解此题）：cos（x）*cosh（x）=-1

!!!using("luopt");
f(x)=cos(x)*cosh(x)+1;
iFind[@f,optmax,2000];  //提高了求解精度

    结果：

-20.42035225104125        -6.825807852273158e-008
-17.27875953208824        1.781003833301043e-008
-14.13716839104647        -9.371081688414051e-011
-10.99554073487547        7.340683616519073e-012
-7.854757438237613        -8.704148513061227e-014
-4.694091132974175        2.364775042451583e-014
-1.875104068711961        2.220446049250313e-016
1.875104068711961         2.220446049250313e-016
4.694091132974175         2.364775042451583e-014
7.854757438237613         -8.704148513061227e-014
10.99554073487547         7.340683616519073e-012
14.13716839104647         -9.371081688414051e-011
17.27875953208824         1.781003833301043e-008
20.42035225104125         -6.825807852273158e-008

14

    （4）luopt::Find(f, luopt::optmode,mode, luopt::optrange,x1min,x1max,x2min,x2max,...,xnmin,xnmax, luopt::opteps,eps)：求方程组的全部解

    f：自定义2*n元函数句柄，不可缺省。格式如下：

f(x1,x2,...,xn,y1,y2,...,yn)=
{
y1= ...,
y2= ...,
... ...,
yn= ...
};

    luopt::optmode,mode：工作模式，取0、1、2、3、... ...。通常，工作模式取值越大，搜到的解越多，但耗时越长。若缺省该参数，工作模式取0。
luopt::optrange,x1min,x1max,x2min,x2max,...,xnmin,xnmax：指定求解区间。若缺省该参数，则所有变量区间均为[-1e50～1e50]。
luopt::opteps,eps：控制精度要求，eps>0。可以缺省该参数，缺省值为1e-6。满足sqrt[(y1*y1+y2*y2+...+yn*yn)/n]<eps。

    返回值：解的个数。

    说明：该函数的解是不稳定的，需要多次运行该函数以获得需要的解。

    例子1：解方程组：

(x-y)^2-3*(x-y) = 10
x^2+2*x*y+y^2 = 9

    代码（最后一个值为误差，下同）：

f(x,y,y1,y2)=
{
y1=(x-y)^2-3*(x-y)-10,
y2=x^2+2*x*y+y^2-9
};
luopt::Find[@f];

    结果：

4.                -1.                0.

0.5               2.5                0.

-2.5              -0.5               0.

1.000000000224727 -4.000000000225374 2.227885009675654e-009

4

    例子2：解方程组：a取-10～10，t取-7～7

-b*sin(a+6*t)+n-40.4945=0
-b*sin(a+7*t)+n-40.5696=0
-b*sin(a+8*t)+n-41.0443=0
-b*sin(a+9*t)+n-41.4190=0

    代码：

!!!using["luopt"];
f(a,b,n,t,y1,y2,y3,y4)=
{
y1=-b*sin(a+6*t)+n-40.4945,
y2=-b*sin(a+7*t)+n-40.5696,
y3=-b*sin(a+8*t)+n-41.0443,
y4=-b*sin(a+9*t)+n-41.4190
};
Find[@f, optmode,5, optrange,-10.,10.,-1e50,1e50,-1e50,1e50,-7.,7.];

    结果（没有验证过该方程组有多少组解，下面是其中一次求解结果）：

8.423278510740619  0.491530082706187   40.94928398718974 -1.077226214994066 5.024295867788081e-015

5.281685857151022  -0.4915300827061788 40.94928398718975 -1.077226214994089 1.123466709944544e-014

4.143092103608209  -0.4915300827069298 40.94928398718954 1.077226214995628  8.519922959291649e-013

1.001499450021447  0.4915300827000679  40.94928398718732 -5.205959092183973 4.316661155959305e-012

-4.143092103780588 0.4915300827292421  40.94928398721459 -1.077226214971404 2.308770421222309e-011

-5.281685858383328 0.4915300826310666  40.949283987277   1.077226215139597  1.209674914149958e-010

-1.001499513894245 -0.4915300850351296 40.94928398852594 5.205959099724587  3.485087466848284e-009

-2.140093572279421 -0.4915300631949757 40.94928396791384 1.077226272463665  1.373241653986354e-008

-1.001499452411673 -0.491530028262163  40.94928398098966 -1.077226219100778 3.989970963218456e-008

9

    （5）luopt::Opt(f, ... ...)：求无约束或约束条件下的n维极小值，可用于解方程（组）

    有些方程（组）用常规方法求解比较困难，须借助于优化函数Opt求解。

    很少有软件可求解以下方程组：

a*exp(-b*7.85)+c*exp(-d*7.85)=28.9
-a/b*exp(-b*7.85)-c/d*exp(-d*7.85)=500
a/b^2*exp(-b*7.85)+c/d^2*exp(-d*7.85)=233
a*exp(-b*1.85)+c*exp(-d*1.85)=20.9

    代码：

!!!using["luopt"];
f(a,b,c,d:f1,f2,f3,f4)=
f1=a*exp(-b*7.85)+c*exp(-d*7.85)-28.9,
f2=-a/b*exp(-b*7.85)-c/d*exp(-d*7.85)-500,
f3=a/b^2*exp(-b*7.85)+c/d^2*exp(-d*7.85)-233,
f4=a*exp(-b*1.85)+c*exp(-d*1.85)-20.9,
sqrt[(f1*f1+f2*f2+f3*f3+f4*f4)/4];
Opt[@f, optwaysimdeep, optwayconfra];

    多次运行，可得以下2组解（a,b,c,d，误差）：

19.08176873960328   -5.36612020972517e-002  -0.1719391618522963 -4.244947571876434e-003 1.70094863414172e-013

-0.1719391618522608 -4.244947571876049e-003 19.08176873960331   -5.366120209725268e-002 7.53644380168212e-015