Lu求解含积分的复杂非线性方程（组）

欢迎访问 Lu程序设计

Lu求解含积分的复杂非线性方程（组）

例子1 解含参变量多重积分的方程组：

解法1：只求一个解

!!!using["math"];
init(::p,q,m,C1,C2,C3,C4,k,g,T3,T2)= p=0.020,q=0.219,m=10369.6,C1=800,C2=2,C3=6,C4=8,k=3,g=4,T3=8.0,T2=12.0;
t_T2(u:a:p,q,m)= a=exp[-(p+q)*u],m*p*(p+q)^2*a/(p+q*a)^2;
t2_T2(t::T2,p,q,g)= gsl_qng[@t_T2,t,T2]/[g*(p+q)];
f(t2,T1,y1,y2:a,b:p,q,m,C1,C2,C3,C4,k,g,T3,T2)= //函数定义
{
    a=exp[-(p+q)*T3], b=m*p*(p+q)^2*a/(p+q*a)^2,
    y1=-C2*gsl_qng[@t2_T2,t2,T2]+C4*b*[1-k*(p+q)] + C3*{k*(p+q)*gsl_qng[@t_T2,T1,T3]+b*k*(p+q)*t2-k*(p+q)*T3},
    y2=C2*T1*T1/[2*g*(p+q)]-k*(p+q)*C3*(t2-T1)-C4*[1-k*(p+q)]
};
main(:x,i,t2,T1)=
{
    o{x=ra1[1,1]}, //修改初值为20,10，还可得到一组解
    i=netn[@f,x],
    o{"\r\n实际迭代次数=",i," t2=",x(0)," T1=",x(1),"\r\n"}
};

结果：

实际迭代次数=4, t2=5.6891908532677187, T1=3.4017560817753707

解法2：全局搜索到4个解（运行时间较长）

!!!using["luopt","math"];
init(::p,q,m,C1,C2,C3,C4,k,g,T3,T2)= p=0.020,q=0.219,m=10369.6,C1=800,C2=2,C3=6,C4=8,k=3,g=4,T3=8.0,T2=12.0;
t_T2(u:a:p,q,m)= a=exp[-(p+q)*u],m*p*(p+q)^2*a/(p+q*a)^2;
t2_T2(t::T2,p,q,g)= gsl_qng[@t_T2,t,T2]/[g*(p+q)];
f(t2,T1,y1,y2:a,b:p,q,m,C1,C2,C3,C4,k,g,T3,T2)=
{
    a=exp[-(p+q)*T3], b=m*p*(p+q)^2*a/(p+q*a)^2,
    y1=-C2*gsl_qng[@t2_T2,t2,T2]+C4*b*[1-k*(p+q)] + C3*{k*(p+q)*gsl_qng[@t_T2,T1,T3]+b*k*(p+q)*t2-k*(p+q)*T3},
    y2=C2*T1*T1/[2*g*(p+q)]-k*(p+q)*C3*(t2-T1)-C4*[1-k*(p+q)]
};
Find[@f];

结果（每行前面的数是解，最后一个数是误差，下同）：

5.689190853227517 3.401756081214042   2.744521617779239e-013
5.03232442248394   -7.261112001867591 2.985639533910636e-012
24.37808017711167 8.27093826714037    3.74498453245254e-012
27.67021280221615 -13.01959458323637 1.01804795021056e-012
145373.6775289801 -775.2875381451365 4.636153342260031e-007
5

例子2 方程组如图所示：

解法1：只求一个解

!!!using["math"];
pp(x::p)=exp{-[(x/p)^2]};
f(_p,q,y1,y2::hpp,p)= //函数定义
{
    p=_p,
    y1=q*gsl_qng[@pp,0.0,p-q]-1.99,
    y2=q*gsl_qng[@pp,0.0,p+q]-2.87
};
main(:x,i,p,q)=
{
    x=ra1[1,1], //1,1是初值
    i=netn[@f,x],
    o{"\r\np=",x(0)," q=",x(1), " 实际迭代次数=",i, "\r\n"}
};

结果：

p=3.201863972597816 q=1.0747323894812955 实际迭代次数=4.

解法2：全局搜索到2个解

!!!using["luopt","math"];
pp(x::p)=exp{-[(x/p)^2]};
f(_p,q,y1,y2::p)=
{
    p=_p,
    y1=q*gsl_qng[@pp,0.0,p-q]-1.99,
    y2=q*gsl_qng[@pp,0.0,p+q]-2.87
};
Find[@f];

结果

3.20186397420115 1.074732389098163 0.
-3.20186397420115 -1.074732389098163 0.
2.

例子3 方程组如图所示：

这个方程组在求解时，可将q消去转换为一元方程求解。但这里不这样做，而是解二元非线性方程组。

全局搜索的Lu代码：

!!!using["luopt","math"];
pp(x::p)=exp{-[(x/p)^2]};
f(_p,q,y1,y2::p)=
{
    p=_p,
    y1=q*gsl_qng[@pp,0.0,0.01]-1.99,
    y2=q*gsl_qng[@pp,0.0,0.015]-2.87
};
Find[@f];

结果：

3.187569597918258e-002 205.5488044354935 0.
-3.187569597918265e-002 205.5488044354935 3.510833468576701e-016
2.

例子4 求解含积分的方程：

取Ф0=45.03，c=0.2277，z1=2100, z2=2730, z3=0；求z4？

Lu代码：

!!!using["luopt","math"];
mvar:
Ф0=45.03, c=0.2277, z1=2100.0, z2=2730.0, z3=0.0;
fun(x)=1-Ф0*exp(-c*x);
f(z4)=gsl_qng[@fun,z3,z4]-gsl_qng[@fun,z1,z2];
iFind[@f];

结果：

-6.320928620456309 0.
827.7602108036889 -4.547473508864641e-013
2

例子5 解方程组：

52.0933-(x1+x2*(1-exp(-(x3*0^x4)))) = 0;
52.2203-(x1+x2*(1-exp(-(x3*9^x4)))) = 0;
53.0575-(x1+x2*(1-exp(-(x3*30^x4)))) = 0;
59.8562-(x1+x2*(1-exp(-(x3*60^x4)))) = 0;

解法1：用优化方法求解

!!!using["luopt","math"];
f(x1,x2,x3,x4)=    //函数定义
{
  [52.0933-(x1+x2*(1-exp(-(x3*0^x4))))]^2+
  [52.2203-(x1+x2*(1-exp(-(x3*9^x4))))]^2+
  [53.0575-(x1+x2*(1-exp(-(x3*30^x4))))]^2+
  [59.8562-(x1+x2*(1-exp(-(x3*60^x4))))]^2
};
ff(x1,x2,x3,x4,y1,y2,y3,y4)=  //函数定义
{
  y1=52.0933-(x1+x2*(1-exp(-(x3*0^x4)))),
  y2=52.2203-(x1+x2*(1-exp(-(x3*9^x4)))),
  y3=53.0575-(x1+x2*(1-exp(-(x3*30^x4)))),
  y4=59.8562-(x1+x2*(1-exp(-(x3*60^x4))))
};
main(:f,x,x1,x2,x3,x4)=
{
  x1=-1.0, x2=-1.0, x3=-1.0, x4=1.0, //初值
  f=Opt[@f, optstarter : &x1,&x2,&x3,&x4,0],  //优化，获得更优的初值
  x=ra1[x1,x2,x3,x4],
  o{"实际迭代次数=",netn[@ff,x]," 方程的解：",x} //拟牛顿法解方程组
};

结果：

52.0933 -0.153590166412464 -6.843318743773691e-002 0.9900707920861045 0.
实际迭代次数=0 方程的解：

52.0933 -0.15359 -6.84332e-002 0.990071

解法2：用Find函数搜索求解

!!!using["luopt"];
f(x1,x2,x3,x4,y1,y2,y3,y4)=
{
    y1=52.0933-(x1+x2*(1-exp(-(x3*0^x4)))),
    y2=52.2203-(x1+x2*(1-exp(-(x3*9^x4)))),
    y3=53.0575-(x1+x2*(1-exp(-(x3*30^x4)))),
    y4=59.8562-(x1+x2*(1-exp(-(x3*60^x4))))
};
Find[@f, optmode,20, optdeep,20, optmax,2000, optrange,-100.0,100.0,-100.0,100.0,-100.0,100.0,-100.0,100.0];

此例求解有一定难度，需多次运行以上代码求解：

52.0933 -0.1346341915383007 -8.137821339044787e-002 0.9556399446745441 3.552713678800501e-015

例子6 解方程组：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

Lu代码：

!!!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, optdeep,5, optmax,1000, optrange,-1e50,1e50,-1e50,1e50,-1e50,1e50,-7.0,7.0];

此例有无穷解，一种可能的结果如下：

-4.143092103627382 0.4915300827138497 40.94928398719285 -1.077226214992139 5.407296744572597e-012
-186.3554660112772 0.4915300826946679 40.94928398716758 -1.07722621506535 1.977372678137429e-011
-8148.289843965879 0.4915300821494349 40.94928398660991 -5.205959091456436 4.018655510508915e-010
-2.140093202816703 -0.491530083327599 40.94928398743284 1.077226214836548 4.382261299437496e-010
-8634741.805115784 0.4915300827124717 40.94928398720907 -5.205959092280351 5.638413030112714e-010
14.7064638097554 0.4915300828176741 40.94928398686785 5.205959093313839 5.790135041881035e-010
-6081.121877884933 0.491530081965423 40.94928398682594 1.077226213077712 8.754386986222557e-010
-109998208.7850408 0.49153008241012 40.94928398776379 -1.077226211365594 8.809486841925538e-010
-44344.58180470909 0.4915300819314578 40.94928398735201 -1.077226234363118 6.031776937013173e-009
9.

例子7 求实数方程复数域内的全部解（所有解）：x^4+2*x*x+10*x-20=0;

本例若用iFind求解，只能获得实数解：

!!!using["luopt"];
f(x)=x^4+2*x*x+10*x-20;
iFind[@f];

结果：

-2.387183302169079 0.
1.331341259641246 0.

用Find求解方程组，可获得复数解（与实数解比较，获得复数解）：

!!!using["luopt","math"];
cf(x,y)= y=x^4+2*x*x+10*x-20;
cc(x,y,y1,y2 : yy)= cf(x$y,&yy), toreal[yy,&y1,&y2];
Find[@cc, optmode,10];

结果：

-0.9978011333203318 2.419443036554579e-009 0.
-2.387183302169079 2.722776498031435e-033 0.
1.331341259641246 0. 0.
0. -2.114742526881128 4.278897697195266e-015
0. 2.114742526881129 1.062884161748515e-014
5

例子8 求解复数方程组：

(2+5i)*x1-x2^(2-3i)-exp(-x1)=0
-(x1^3)+x1*x2-exp(-x2)=0

Lu代码：

!!!using["luopt","math"];
cf(x1,x2,y1,y2)=
{
y1=(2+5i)*x1-x2^(2-3i)-exp(-x1),
y2=-(x1^3)+x1*x2-exp(-x2)
};
cc(x11,x12,x21,x22,y11,y12,y21,y22 : y1,y2)= cf(x11$x12,x21$x22,&y1,&y2), toreal[y1,&y11,&y12], toreal[y2,&y21,&y22];
Find[@cc, optmode,10];

此例有无穷解，一种可能的结果如下：

-2.866006054064824e-002     -1.195892184644717e-002     1.330587481554971     -8.406178015402372     1.387778780781446e-017
8.343862612043082e-002     -0.1745973157291573     0.3407059687267944     -3.686653995031711     8.441528768080324e-017
8.66377312425692e-003     5.927370162532637e-002     2.009821712467028     -0.9744617637089106     1.666779512946625e-016
0.1575094449966265     -6.233406303054641e-003     -0.542816052357006     -10.90400468060004     9.601098942365283e-016
0.6041474952500348     0.8834204997637829     -2.907695682419187e-002     0.1865018284480821     1.221995758683528e-014
0.3504034061227544     -0.2581172046401701     0.9031492305415145     0.2062068702236728     3.335961918797048e-014
-1.566817874139894     -0.1557788051049168     2.383001997667281     0.5222697207582948     1.165852988550321e-013
0.3248221986169174     -0.8164562030908124     8.588879599544659e-002     0.2894415311814788     8.437690236792216e-011
0.2558786571311035     -0.6418806813542586     -2.724268338847976     -22.2624989134893     7.712800937949101e-010
0.3678174169101899     -0.2447841797871731     -1.987462607248696     -16.56633273776007     6.305964240795578e-009
10
　

例子9 求解方程组：

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

此题看似简单，但一般的解方程算法恐怕都很难求解，以下是Lu的优化解法：

!!!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];
Opt1[@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

Lu求解含积分的复杂非线性方程（组）相关推荐

非线性方程(组)：高维方程解法
非线性方程的高维情形和一维情形既有相似处也有差异.首当其中的区别即在高维情形中不再存在介值定理,从而使得二分法不再可推广到高维.不过,仍然有许多方法可以推广. 1. 不动点迭代(高维) 寻找方程 $\ ...
L1正则化降噪，对偶函数的构造，求解含L1正则项的优化问题，梯度投影法
L1正则化降噪,对偶函数的构造,求解含L1正则项的优化问题,梯度投影法本文主要实现L1正则化降噪,L2 正则化降噪的文章在: https://blog.csdn.net/IYXUAN/article ...
matlab采用粒子群优化算法求解含压缩储能设备的综合能源系统运行优化
matlab采用粒子群优化算法求解含压缩储能设备的综合能源系统运行优化. 结果包含储能设备24时出力,内燃机发电和发热出力,电制冷机出力等. 代码包含相关注释,方便对算法进行改进. 附相关参考文献. ...
（程序）MALTAB求解含未知数的矩阵逆
MATLAB可以求解含未知数的矩阵的逆,下面用一个例子进行说明: 例:对于下面这样的矩阵A,要求它的逆 MATLAB程序 syms s w0 kp kd % 定义未知量 A = [s+3*w0,-1, ...
非线性方程(组)：计算基本理论
1. 非线性方程(组)及其解对于非线性函数 $f(x)$ ,方程 $f(x)=0$ 为非线性方程( $f:\mathbb{R} \rightarrow \mathbb{R}$ ). 对于n维向量 $ ...
非线性方程(组)：一维非线性方程（一）二分法、不动点迭代、牛顿法 [MATLAB]...
1. 二分法(Bisection) 1) 原理 [介值定理] 对于连续的一元非线性函数,若其在两个点的取值异号,则在两点间必定存在零点. [迭代流程] 若左右两端取值不同,则取其中点,求其函数值,取中 ...
matlab求解含两个累加公式的方程,使用fsolve 函数求解含两个参数的多个方程问题，希望有会的人帮...
该楼层疑似违规已被系统折叠隐藏此楼查看此楼使用fsolve 函数求解四个方程,a,b是参数,t(1)_t(4)为变量,求解t(1)随a,b变化 function M= SHI_2( t ) a = ...
【无人机】基于matlab蚁群算法求解含危险源的无人机路径规划【含Matlab源码 2059期】
一.无人机简介 0 引言随着现代技术的发展,飞行器种类不断变多,应用也日趋专一化.完善化,如专门用作植保的大疆PS-X625无人机,用作街景拍摄与监控巡察的宝鸡行翼航空科技的X8无人机,以及用作水下 ...
maple 假设_一道含着假设的不等式组Maple解法
下面的问题具有一定的普遍性,带着假设的方程往往求解会出现空解或者是漏解警告,不能提供更多的信息,大佬acer 提供了一个有益的尝试的解法,主要是利用化简命令的参数化选项进行求解,很值得去体会. 问题描 ...

Lu求解含积分的复杂非线性方程（组）

Lu求解含积分的复杂非线性方程（组）相关推荐

最新文章

热门文章