當(dāng)我們面對(duì)一個(gè)非線性的ODE時(shí)，我們要怎么解決它？

推薦看這篇介紹文章：http://www.phys.uconn.edu/~rozman/Courses/P2400_15S/downloads/averaging.pdf

這里寫下我的理解：

以一個(gè)單自由度的mass-damper-spring系統(tǒng)為例，但是這里的阻尼不是關(guān)于速度線性的，而是關(guān)于速度的三次非線性的，其控制方程和初始條件如下所示：

1. 首先用數(shù)值方法得到數(shù)值解作為對(duì)比，將二階的ODE寫成兩個(gè)一階的ODEs：

該ODEs可以直接用Matlab的ode45求解，代碼如下：

% numerical solution by ode45
epsilon = 0.2;
tspan = 0:0.01:50;
y0 = [1; 0];
[~, y] = ode45(@(t,y) equations(y, epsilon), tspan, y0);
%  each row in y corresponds to the solution at the value returned
%  in the corresponding row of t.
plot(tspan, y(:,1),'linewidth',1.5)
grid on

function dydt = equations(y, epsilon) 
% oscillator with nonlinear friction
% y(1) ---> x
% y(2) ---> dx/dt
y = y(:);
dydt = zeros(size(y));
dydt(1) = y(2);
dydt(2) = -epsilon*y(2)^3-y(1);
end

　結(jié)果如下圖：

　

表示位移x(t)隨時(shí)間的變化

2.用攝動(dòng)法（perturbation）求近似解析解

首先將x(t)按epsilon的冪次展開(kāi)

（為什么要這樣展開(kāi)？我的理解：把解分成幾部分之和，每一部分的幅值大小由epsilon刻畫，因?yàn)閑psilon是一個(gè)小量，顯然x0(t)對(duì)x(t)的貢獻(xiàn)最主要，x1(t)次之，貢獻(xiàn)依次遞減。也就是說(shuō)認(rèn)為x(t)主要是按x(t)變化，受到了一些攝動(dòng)x1(t)，攝動(dòng)x1(t)對(duì)x(t)的影響大小與小量epsilon同階。這種展開(kāi)應(yīng)該主要應(yīng)用在弱非線性的系統(tǒng)，即認(rèn)為解x(t)仍然主要按對(duì)應(yīng)線性系統(tǒng)的解x0(t)變化，再加上一點(diǎn)非線性的影響（攝動(dòng)）x1(t),x2(t)，這些攝動(dòng)的影響是小量的。因?yàn)槭菑男×縠psilon對(duì)解的影響來(lái)看待問(wèn)題的，所以這種方法主要適用于弱非線性系統(tǒng)）

按這種方式展開(kāi)之后，為了與原初始條件協(xié)調(diào)一致，這些分量的初始條件應(yīng)該為

將展開(kāi)式（5）式代入原方程（1）式，并讓同一epsilon的冪次的合并，可以得到一系列線性的方程，注意看到由于非線性（這里是三次方）的存在，展開(kāi)的時(shí)候顯然很麻煩，這也是攝動(dòng)法的一個(gè)缺點(diǎn)

這里只寫到epsilon的一次項(xiàng)，一般來(lái)說(shuō)對(duì)于近似解這就足夠的，而且更高次的寫出來(lái)顯然很麻煩

注意到第一個(gè)方程，x0，是線性的，其對(duì)x的貢獻(xiàn)是epsilon的零次的，也就是最主要的，接下來(lái)epsilon的一次冪的解x1對(duì)x的貢獻(xiàn)與epsilon的一次同階，表示攝動(dòng)。在下面的分析中我們這考慮x1，其他攝動(dòng)量都是很小的（其大小程度由epsilon的冪次衡量）

線性部分（主要部分）的滿足初始條件的解是

對(duì)于epsilon的一階項(xiàng)，注意到x0已經(jīng)求出來(lái)了，所以是線性的方程。也就是說(shuō)用攝動(dòng)法我們可以得到對(duì)每一個(gè)epsilon的冪次的線性方程，這是攝動(dòng)法的優(yōu)點(diǎn)，把非線性方程分成幾個(gè)線性的方程求解。但是要注意到等式右邊一般來(lái)說(shuō)是關(guān)于時(shí)間t的復(fù)雜的表達(dá)式（一般是三角函數(shù)），所以即使是要求解一個(gè)線性微分方程的解析解也是麻煩的

以（7）式子為例，

至此我們可以得到一個(gè)近似的解析解：

將該結(jié)果與上面用ode45得到的數(shù)值解作對(duì)比

代碼如下

% numerical solution by ode45
epsilon = 0.2;
tspan = 0:0.01:50;
y0 = [1; 0];
[~, y] = ode45(@(t,y) equations(y, epsilon), tspan, y0);
%  each row in y corresponds to the solution at the value returned
%  in the corresponding row of t.
plot(tspan, y(:,1),'b','linewidth',1.5)
grid on
hold on
% perturbation solution
x1 = cos(tspan)+epsilon*(3/32*sin(tspan)-1/32*sin(3*tspan)-3/8*tspan.*sin(tspan));
plot(tspan,x1, 'r', 'linewidth', 1.5)

legend('ode45(benchmark)','perturbation')



function dydt = equations(y, epsilon) 
% oscillator with nonlinear friction
% y(1) ---> x
% y(2) ---> dx/dt
y = y(:);
dydt = zeros(size(y));
dydt(1) = y(2);
dydt(2) = -epsilon*y(2)^3-y(1);
end

　　

可見(jiàn)攝動(dòng)法的解只有在開(kāi)始的一段很短的時(shí)間是對(duì)的，很快后面就不對(duì)了，會(huì)發(fā)散（變得很大）。所以攝動(dòng)法是有問(wèn)題的。問(wèn)題出現(xiàn)在什么地方？讓我們看一下用攝動(dòng)法得到的解

顯然是由于這一項(xiàng)的存在才會(huì)使得解會(huì)發(fā)散（grow undounded as t increases），這一項(xiàng)叫做長(zhǎng)期項(xiàng)（secular term）。雖然一般翻譯成長(zhǎng)期項(xiàng)，但是我想英文想表達(dá)的是（not bound by monastic vows or rules）的意思，（翻譯成病態(tài)項(xiàng)或許更好？），而不是長(zhǎng)期存在，因?yàn)閯e的項(xiàng)也長(zhǎng)期存在（比如sin(t), sin(3t)），所以接下里不用長(zhǎng)期項(xiàng)這種翻譯而是直接用secular term

secular term之所以會(huì)出現(xiàn)是因?yàn)樵冢?1）式的右邊（激勵(lì)）有這一部分，注意到其激勵(lì)頻率等于左邊系統(tǒng)的固有頻率，即發(fā)生了共振。

這是攝動(dòng)法的一個(gè)常見(jiàn)缺點(diǎn)（或者說(shuō)致命缺點(diǎn)？），至于為什么會(huì)出現(xiàn)，什么情況下可以避免，這里暫時(shí)不深入。總之只需要知道攝動(dòng)法求解時(shí)有時(shí)候會(huì)由于（人為地引入）共振的出現(xiàn)，導(dǎo)致解有secular term。

3. 平均法

平均法可以用來(lái)求解下面這一類關(guān)于一階導(dǎo)數(shù)非線性的非線性O(shè)DE：

上一節(jié)用攝動(dòng)法考慮的例子是一個(gè)特例：

平均法考慮（13）式具有如下形式的解和速度解，注意不僅假設(shè)了位移的形式，同時(shí)也假設(shè)了速度具有和線性時(shí)相同的形式：

為什么要這么考慮？因?yàn)楫?dāng)非線性消失時(shí)，退化的線性O(shè)DE的解具有（15）式和（16）式的形式，不過(guò)此時(shí)的幅值和相位是由初始條件決定的常數(shù)；當(dāng)系統(tǒng)有了一點(diǎn)弱非線性時(shí)，我們可以樂(lè)觀地假設(shè)方程的解仍然具有（15）式的形式，不過(guò)此時(shí)幅值和相位是會(huì)隨時(shí)間變化的函數(shù)，并且認(rèn)為幅值和相位是緩慢變化的。為什么認(rèn)為是緩慢變化的？因?yàn)榫徛兓脑挿匠痰慕獠糯笾碌鼐哂泻途€性時(shí)候的解一樣的形式，只不過(guò)由于非線性帶來(lái)幅值和相位的緩慢變化。注意這里頻率仍然是固定的。若非線性出現(xiàn)在會(huì)影響頻率的剛度項(xiàng)，又該怎么做呢？

這是平均法的核心，或者說(shuō)亮點(diǎn)吧，就是洞見(jiàn)問(wèn)題的解具有（15）和（16）式的形式。

將（15）式求導(dǎo)代入（16）式得到一個(gè)方程，將（16）式求導(dǎo)代入（1）式，得到兩個(gè)耦合的ODEs

整理成

到這一步為止都是精確的，不過(guò)是假設(shè)了解的形式，但沒(méi)有作任何的近似或者截?cái)唷?/p>

接下來(lái)為了求解（20）和（21）式，要做一點(diǎn)假設(shè)：

因?yàn)槭侨醴蔷€性，即epsilon是個(gè)小量，所以認(rèn)為幅值和相位的變化是緩慢的，即幅值和相位的導(dǎo)數(shù)是小量。這樣的話，在振動(dòng)的一個(gè)周期內(nèi)（這里頻率是1，所以周期是2pi），（20）式和（21）式的方程右邊幾乎保持不變是一個(gè)常數(shù)，這樣就可以用它們?cè)谝粋€(gè)周期內(nèi)的平均作近似。（這樣的近似還是不是很理解）

（20）和（21）式變?yōu)橄旅娴男问剑瞧骄姆匠?/p>

注意到積分后時(shí)間變量就消去了（消除快變量），（不是很理解）

對(duì)（23）式和（24）式右邊的積分計(jì)算，需要用到三角恒等式

最終（23）和（24）式變?yōu)?/p>

（26）式和（27）式可以較輕易地求解：

最終解為下面的形式，其中a(0)由初始條件確定為x(0)

和數(shù)值解作比較發(fā)現(xiàn)是較為精確地！而且和數(shù)值解相比具有解析形式的優(yōu)勢(shì)。

代碼如下：

% numerical solution by ode45
epsilon = 0.2;
tspan = 0:0.01:50;
y0 = [1; 0];
[~, y] = ode45(@(t,y) equations(y, epsilon), tspan, y0);
%  each row in y corresponds to the solution at the value returned
%  in the corresponding row of t.
plot(tspan, y(:,1),'b','linewidth',2)
grid on
hold on
% perturbation solution
x1 = cos(tspan)+epsilon*(3/32*sin(tspan)-1/32*sin(3*tspan)-3/8*tspan.*sin(tspan));
plot(tspan,x1, 'r', 'linewidth', 1.5)
% averaged solution
x2 = cos(tspan)./sqrt(3/4*epsilon*tspan+1/1^2);
plot(tspan,x2,'y--')
legend('ode45(benchmark)','perturbation','averaged')



function dydt = equations(y, epsilon) 
% oscillator with nonlinear friction
% y(1) ---> x
% y(2) ---> dx/dt
y = y(:);
dydt = zeros(size(y));
dydt(1) = y(2);
dydt(2) = -epsilon*y(2)^3-y(1);
end

　　

4. 接下來(lái)考慮另一個(gè)系統(tǒng)，這是一個(gè)著名的非線性系統(tǒng)，叫做van der pol振子，具有（13）式的形式，因此也可以用平均法求近似解析解

求數(shù)值解（ode45）時(shí)仍然要先轉(zhuǎn)化為一階的ODEs

求近似解析解時(shí)，和上面一樣的思路和步驟，首先假設(shè)解的形式為：

第一式求導(dǎo)代入第二式，第二式求導(dǎo)代入方程，得到

整理：

接下來(lái)做平均：

計(jì)算積分同樣要用到三角恒等式

最后得到

求解ODE（實(shí)際還很難求解的）：

和數(shù)值解作比較，在前面一段時(shí)間有點(diǎn)誤差，后面對(duì)得還好：

代碼如下

% numerical solution by ode45
epsilon = 0.1;
tspan = 0:0.01:100;
y0 = [1; 0];
[~, y] = ode45(@(t,y) equations(y, epsilon), tspan, y0);
%  each row in y corresponds to the solution at the value returned
%  in the corresponding row of t.
plot(tspan, y(:,1),'b','linewidth',2)
grid on
hold on
% averaged solution
x2 = (2-exp(-epsilon*tspan)).*cos(tspan);
plot(tspan,x2,'r','linewidth',2)
legend('ode45(benchmark)','averaged')



function dydt = equations(y, epsilon) 
% oscillator with nonlinear friction
% y(1) ---> x
% y(2) ---> dx/dt
y = y(:);
dydt = zeros(size(y));
dydt(1) = y(2);
dydt(2) = -epsilon*(y(1)^2-1)*y(2)-y(1);
end

　　

總結(jié)：

平均法的兩個(gè)步驟：

一是假設(shè)解的形式，這個(gè)形式由非線性激勵(lì)和/或非線性阻尼中的小參數(shù)取零時(shí)確定，即系統(tǒng)在無(wú)阻尼無(wú)激勵(lì)時(shí)的精確解形式，然后令幅值和相位變?yōu)闀r(shí)間的函數(shù)

二是認(rèn)為幅值和相位的變化是小量，在一個(gè)周期內(nèi)不變，所以可以對(duì)時(shí)間在一個(gè)周期上上作平均，去掉快變量

練習(xí)：
1.Duffing

我的解答：

代碼如下

% numerical solution by ode45
epsilon = 0.2;
tspan = 0:0.01:25;
y0 = [1; 0];
[~, y] = ode45(@(t,y) equations(y, epsilon), tspan, y0);
%  each row in y corresponds to the solution at the value returned
%  in the corresponding row of t.
plot(tspan, y(:,1),'b','linewidth',2)
grid on
hold on
% averaged solution
x2 = cos((1+3/8*epsilon)*tspan);
plot(tspan,x2,'r','linewidth',2)
legend('ode45(benchmark)','averaged')



function dydt = equations(y, epsilon) 
% oscillator with nonlinear friction
% y(1) ---> x
% y(2) ---> dx/dt
y = y(:);
dydt = zeros(size(y));
dydt(1) = y(2);
dydt(2) = -epsilon*y(1)^3-y(1);
end

3.

我的解答

數(shù)值解對(duì)比

　代碼：

% numerical solution by ode45
epsilon = 0.03;
tspan = 0:0.01:25;
y0 = [1; 0];
[~, y] = ode45(@(t,y) equations(y, epsilon), tspan, y0);
%  each row in y corresponds to the solution at the value returned
%  in the corresponding row of t.
plot(tspan, y(:,1),'b','linewidth',3)
grid on
hold on
% averaged solution
x2 = sqrt(2)./(4+5*epsilon*tspan).^(1/4).*cos(tspan);
plot(tspan,x2,'r','linewidth',1.5)
legend('ode45(benchmark)','averaged')



function dydt = equations(y, epsilon) 
% oscillator with nonlinear friction
% y(1) ---> x
% y(2) ---> dx/dt
y = y(:);
dydt = zeros(size(y));
dydt(1) = y(2);
dydt(2) = -epsilon*y(2)^5-y(1);
end

　2.

看起來(lái)有點(diǎn)麻煩，罷了

后續(xù)應(yīng)對(duì)平均法作深入了解，是求解非線性微分方程的近似解析解的主要方法

　

總結(jié)

以上是生活随笔為你收集整理的非线性微分方程的平均法的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。