应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 专业学科区 » 数学 » 计算数学 » 请教一个非线性方程非线性项转化为多项式的方法
321/4返回列表1234下一页
查看: 3000  |  回复: 31
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

businiao412

金虫 (正式写手)

[求助] 请教一个非线性方程非线性项转化为多项式的方法 已有2人参与

RT,楼主近日做振动问题时遇到一个二阶非线性微分方程如图片1，看文献中其有如下解法:在运用摄动法之前，先将其中的非线性项（1-q）3/2运用三阶最小二乘拟合的方法在区间1-Q<q<1上进行展开如图片2，及图片3，请问这是怎么做到的呢？我查了很多资料就是找不到这个方法的出处，本人数学知识较为薄弱，望不吝赐教。
图片1-3如下图所示：
请教一个非线性方程非线性项转化为多项式的方法
方程1.png


请教一个非线性方程非线性项转化为多项式的方法-1
2.png


请教一个非线性方程非线性项转化为多项式的方法-2
3.png

[ Last edited by businiao412 on 2014-8-10 at 15:16 ]
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼 2014-08-10 15:11:01
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
回帖支持 ( 显示支持度最高的前 50 名 )

guojiashun

铁杆木虫 (正式写手)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与，应助指数 +1
businiao412: 金币+15, ★★★★★最佳答案, 衷心感谢您细致透彻的讲解，有种豁然开朗的感觉，我按照您的思路试做一下，有问题再向您请教。 2014-08-11 14:57:17
对你的问题理解不是十分准确。仅仅就我的理解说一下：
      摄动法的基本思想是把响应q看做多个不同量级的摄动的叠加，即：q(t,ε )=u0(t)+ε*u1(t)+ε ^2*u2(t)+......【1】（多尺度法还要把考虑不同的时间尺度，q(t,ε)=u0(t1,t2,...)+ε*u1(t1,t2,...)+ε ^2*u2(t1,t2,...)+.....），然后把q的近似表达式【1】代入原方程，通过整理和推导，根据方程两边 小参数 ε 的同次幂系数相等的条件，把原方程转化为分别关于u0、u1、u2.....的微分方程，然后分别逐次求得各摄动u0、u1、u2........最后，叠加可得方程q的近似解析解。
      这个非线性微分方程的非线性项是（1-q）^3/2，而在运用摄动法把式【1】代入原方程中进行整理时，由于（1-q）^3/2这一项带有分数次幂，故【1】代入后其展开是非常困难的，即使展开其中也带有 ε 的分数次幂项，也就无法利用 ε 的同次幂系数相等的条件来求解（后面的套路步骤走不下去了）。所以，一般在求解此类非线性微分方程时，用多项式来近似表达诸如（1-q）^3/2这种难于处理的项，这样当式【1】代入后展开、整理就变得很简便了。类似的情况如齿轮动力学理论中，一个主要的非线性因素是齿侧间隙，间隙作用方式用分段函数来表达，而分段函数的存在也给摄动法、多尺度法等的运用带来困难，数学上难于展开处理，故也有学者用多项式来近似这个分段函数（当然还有其它处理方法，这个还要视造成困难的项的具体函数性质）。
      这里提前预估了一个拟合区间q∈[1-Q, 1]，Why？第一，方程首先要有稳态解；第二，也是为了数学处理上的便利性。（1）根据（1-q）^3/2可知，q必须小于等于1，否则出现复数项，方程无稳态解，故已知q的上限；（2）采用多项式来拟合（1-q）^3/2，取多少阶合适？我们当然是希望阶数越低越好。阶数越高，式【1】代入后的数学展开越复杂繁琐。根据数学知识易知拟合区间越小，就可以用越低的阶数来获得足够的拟合精度，所以根据已知条件提前划定一个尽量小的区间十分有益和必要。接下来要确定q的下限，其中Q是q所能预期到的最大的峰值（上限）和谷值（下限）之间的差值（peak-to-peak value，我想这里的Q可以通过数值求解、实验或者实际观测得到），那么已知q的上限是1，下限自然就是1-Q了。
      最后，文中确定用三次多项式来拟合，接着就要确定多项式的系数a0、a1、a2、a3。这里采用“最小二乘法”来确定，最小二乘法是一种数值方法，必须基于已知的数据进行拟合。上面已经得到了拟合区间q∈[1-Q, 1]，那么在这个区间内取一些点，如1-Q, 1-9*Q/10, 1-8*Q/10, 1-7*Q/10, ......, 1，依次代入（1-q）^3/2中，可得相应值Q^3/2, (9*Q/10)^3/2, (8*Q/10)^3/2, (7*Q/10)^3/2, ......, 0，已知了函数（1-q）^3/2的一组自变量和因变量，套用“最小二乘法”的套路解得多项式的系数a0、a1、a2、a3也就不难了吧？
      总结，不论是摄动法，还是多项式拟合，都是无限逼近思想的体现。摄动法本身得到的就是方程的无限近似解，因此用多项式无限近似方程中复杂项也就变得无可厚非了。

PS：周日下午无事，多敲了些废话，纯属个人见解，限于水平，有误之处请各位指正。
赞一下(2人) 回复此楼
3楼2014-08-10 19:43:10
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

herowolf

木虫 (正式写手)

【答案】应助回帖

感谢参与，应助指数 +1
(1+x)ⁿ = 1 + n·x/1! + n(n-1)x2/2!  + n(n-1)(n-2)x3/3! + . . . + n(n-1)...(n-k+1)xⁿ/n! + . .

是不是这样做的？ -x -> x  ,  n ->3/2

third order 只去前三项。
赞一下(1人) 回复此楼
2楼2014-08-10 16:53:53
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

guojiashun

铁杆木虫 (正式写手)
引用回帖:
4楼: Originally posted by businiao412 at 2014-08-11 14:59:13
十分感谢您细致的解答，对于数学小白的我还有很长的路要走，希望以后不吝赐教。...

我也是小白，只是对这方面稍微熟悉一点，以后多交流，共同进步。
赞一下(1人) 回复此楼
6楼2014-08-11 15:50:17
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

feixiaolin

荣誉版主 (文坛精英)

优秀版主
建议  
尝试 4个控制点的里米兹多项式多项式。
赞一下(1人) 回复此楼
7楼2014-08-11 16:08:40
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

guojiashun

铁杆木虫 (正式写手)
引用回帖:
12楼: Originally posted by businiao412 at 2014-08-13 22:45:37
syms Q
x=;
y=;
P=polyfit(x,y,3)我用这样的表达式，但是一直出错，如下信息：
Error using ==> ones
Trailing string input must be a valid numeric class name.

Error in ==> polyfit at 63
V(: ...

那是因为matlab的ones函数中不支持符号运算，合法的数据类包括： 'double', 'single', 'int8', 'uint8', 'int16', 'uint16', 'int32', 'uint32', 'int64', or 'uint64'，不包括符号类：sym（详见关于ones函数的help文件）。matlab的符号运算功能还是比较弱，你可以试试用maple。

由于我没装maple，我直接用具体数值替你做了几组验证，证明我说的基本是正确的：
（1）令Q=1，则x=[0,0.1,0.2,......,0.9,1]，对应的y=-(1.-x).^(3/2)=[-1,-0.853814968245462,-0.715541752799933,....,0]，p=polyfit(x,y,3)=-0.264788988716828  -0.187032536397084   1.452896005112752  -0.998323698149733，根据论文知：
a3=-8/(33*Q^(3/2))=-0.242424242424242；
a2=-(72*Q-56)/(77*Q^(3/2))=-0.207792207792208；
a1=(24*Q^2+144*Q-56)/(77*Q^(3/2))=1.454545454545455；
a0=(8*Q^3-360*Q^2-1080*Q+280)/(1155*Q^(3/2))=-0.997402597402597；

（2）令Q=0.5，则x=[0.5,0.55,0.6,......,1]，对应的y=-(1.-x).^(3/2)=[-0.3535,-0.3018,-0.2529,......,-0.0111,0]，p=polyfit(x,y,3)=-0.748936358020765   0.858900587453328   0.730154298636324  -0.839145979819353，则
a3=-8/(33*Q^(3/2))=-0.685679302968773；
a2=-(72*Q-56)/(77*Q^(3/2))=0.734656396037972；
a1=(24*Q^2+144*Q-56)/(77*Q^(3/2))=0.808122035641769；
a0=(8*Q^3-360*Q^2-1080*Q+280)/(1155*Q^(3/2))=-0.854650274057507；

结论：可见在区间[1-Q,1]内取点进行拟合得到的多项式系数a3,a2,a1,a0的值与利用论文中a3,a2,a1,a0公式得到的数值是能够对上的。中间有些误差，主要原因一是取点的疏密对拟合精度有影响；二是论文中得到的表达式也是近似的。
赞一下(1人) 回复此楼
13楼2014-08-14 01:43:25
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

guojiashun

铁杆木虫 (正式写手)
内容已删除
赞一下(1人) 回复此楼
15楼2014-08-14 01:55:19
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

guojiashun

铁杆木虫 (正式写手)
引用回帖:
18楼: Originally posted by businiao412 at 2014-08-14 02:40:30
您解答的太给力了，再次表示感谢，您用数据验证了文献的方法的正确性，希望能够用您说的解方程组的形式解出来，不过我心里也有同样的疑问，解出来的结果会不会跟去的数值点的大小跟疏密有关系？...

解方程组的方法应该和点取多少没关系，四个未知量，凑够四个方程就行。我试了一下：取x=[1-Q,1-6*Q/10,1-3*Q/10,1]，y=[-Q^(3/2),-(6*Q/10)^(3/2),-(3*Q/10)^(3/2),0]，
X=[x(1)^3,x(1)^2,x(1),1;
      x(2)^3,x(2)^2,x(2),1;
      x(3)^3,x(3)^2,x(3),1;
      x(4)^3,x(4)^2,x(4),1;]
a=inv(X)*y'
=25/7/Q^3*conj(Q^(3/2))-5/3/Q^3*3^(1/2)*5^(1/2)*conj(Q^(3/2))+10/21/Q^3*3^(1/2)*10^(1/2)*conj(Q^(3/2))
   15/14/Q^3*(-10+3*Q)*conj(Q^(3/2))-1/6/Q^3*(-30+13*Q)*3^(1/2)*5^(1/2)*conj(Q^(3/2))+2/21/Q^3*(-15+8*Q)*3^(1/2)*10^(1/2)*conj(Q^(3/2))
   3/14/Q^3*(50-30*Q+3*Q^2)*conj(Q^(3/2))-1/6/Q^3*(30-26*Q+3*Q^2)*3^(1/2)*5^(1/2)*conj(Q^(3/2))+2/21/Q^3*(15-16*Q+3*Q^2)*3^(1/2)*10^(1/2)*conj(Q^(3/2))
   -1/14/Q^3*(50-45*Q+9*Q^2)*conj(Q^(3/2))+1/6/Q^3*(10-13*Q+3*Q^2)*3^(1/2)*5^(1/2)*conj(Q^(3/2))-2/21/Q^3*(5-8*Q+3*Q^2)*3^(1/2)*10^(1/2)*conj(Q^(3/2))

根据论文知道Q必然大于0，所以conj(Q^(3/2))=Q^(3/2)，这样整理符号运算所得第一项：a3=(25/7-5*sqrt(15)/3+10*sqrt(30)/21)/Q^(3/2)=-0.2753/Q^(3/2)≈  -8/(33*Q^(3/2))，其它不再逐一整理。符号运算得到的a3表达式和论文中还是有误差，这应该是matlab矩阵求逆算法的误差造成的。四元一次方程组应该也不难解，如果手工解是不是应该完全精准？
赞一下(1人) 回复此楼
19楼2014-08-14 14:27:18
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

guojiashun

铁杆木虫 (正式写手)
我用最小二乘法做了一下：
syms Q
x=[1-Q:0.001*Q:1];
y=-(1.-x).^(3/2);
N=length(x);
syms A
for n=1:N
       for m=1:4
           A(m,n)=x(n)^(m-1);
       end
end
a=(A*A')\(A*y');
其中a(1)=a0;  a(2)=a1; a(3)=a2; a(4)=a3;然后得到很长很长很长一串a关于Q的表达式，人工化简基本不可能了.............然后用eval函数化简最简单的 a3 项，得
eval(a(4))=-1230776300378990335736947704895/5070602400912917605986812821504/conj(Q)^3*conj(Q^(3/2))= -0.2427/Q^(3/2)= -8/(33*Q^(3/2))
但是其它几项太复杂了，化简不动.........
赞一下(1人) 回复此楼
29楼2014-08-14 19:03:52
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

guojiashun

铁杆木虫 (正式写手)
把取点密度减小：x=[1-Q:0.1*Q:1]; 可得：
eval(a(4))=-23980501898594795/90564573756653568/conj(Q)^3*conj(Q^(3/2))
             =-0.2648/Q^(3/2);
eval(a(3))=-1/1716*conj(Q^(3/2))*(-2997562737324347/2199023255552-12000*conj(Q)-500*conj(Q)*2^(1/2)+3008*conj(Q)*5^(1/2)+1116*conj(Q)*10^(1/2)+721*conj(Q)*7^(1/2)*10^(1/2)-704*conj(Q)*2^(1/2)*5^(1/2)-381*conj(Q)*3^(1/2)*10^(1/2)+624*conj(Q)*3^(1/2)*5^(1/2))/conj(Q)^3
             =-(1.6841e+003*Q-1.3631e+003)/(1716*Q^(3/2))
             =-(0.9814*Q-0.7943)/Q^(3/2) ≈ -(72*Q-56)/(77*Q^(3/2))

后面实在整不动了，说明最小二乘可定是可行的，但是这个论文作者采用何等手段获得如此整齐精准的表达式呢？
赞一下(1人) 回复此楼
30楼2014-08-14 19:28:03
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
普通回帖

businiao412

金虫 (正式写手)
引用回帖:
3楼: Originally posted by guojiashun at 2014-08-10 19:43:10
对你的问题理解不是十分准确。仅仅就我的理解说一下：
      摄动法的基本思想是把响应q看做多个不同量级的摄动的叠加，即：q(t,ε )=u0(t)+ε*u1(t)+ε ^2*u2(t)+......【1】（多尺度法还要把考虑不同的时间尺度， ...

十分感谢您细致的解答，对于数学小白的我还有很长的路要走，希望以后不吝赐教。
赞一下 回复此楼
4楼2014-08-11 14:59:13
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

businiao412

金虫 (正式写手)
引用回帖:
2楼: Originally posted by herowolf at 2014-08-10 16:53:53
(1+x)ⁿ = 1 + n·x/1! + n(n-1)x2/2!  + n(n-1)(n-2)x3/3! + . . . + n(n-1)...(n-k+1)xⁿ/n! + . .

是不是这样做的？ -x -> x  ,  n ->3/2

third order 只去前三项。

谢谢您热心回复，您说的这个我考虑了，跟文章中的不太一样，再次感谢。
赞一下 回复此楼
5楼2014-08-11 15:14:27
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

businiao412

金虫 (正式写手)
引用回帖:
6楼: Originally posted by guojiashun at 2014-08-11 15:50:17
我也是小白，只是对这方面稍微熟悉一点，以后多交流，共同进步。...

好的，多交流，多向您学习。
赞一下 回复此楼
8楼2014-08-11 16:33:54
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

businiao412

金虫 (正式写手)
引用回帖:
3楼: Originally posted by guojiashun at 2014-08-10 19:43:10
对你的问题理解不是十分准确。仅仅就我的理解说一下：
      摄动法的基本思想是把响应q看做多个不同量级的摄动的叠加，即：q(t,ε )=u0(t)+ε*u1(t)+ε ^2*u2(t)+......【1】（多尺度法还要把考虑不同的时间尺度， ...

再请教您一下，因为不能确定Q的取值，在用matlab拟合过程中，已知参量能否包含Q，写成这样x=[1-Q,1-0.9*Q,1-0.8*Q....],y=[Q^3/2,(0.9)^3/2,(0.8*Q)3/2...]？
赞一下 回复此楼
9楼2014-08-13 11:34:35
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

guojiashun

铁杆木虫 (正式写手)
引用回帖:
9楼: Originally posted by businiao412 at 2014-08-13 11:34:35
再请教您一下，因为不能确定Q的取值，在用matlab拟合过程中，已知参量能否包含Q，写成这样x=,y=？...

可以，matlab有符号运算，或者用maple也可以。
赞一下 回复此楼
10楼2014-08-13 16:23:32
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 businiao412 的主题更新
321/4返回列表1234下一页
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭