當我們面對一個非線性的ODE時，我們要怎么解決它？

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

這里寫下我的理解：

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

1. 首先用數值方法得到數值解作為對比，將二階的ODE寫成兩個一階的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

　結果如下圖：

　

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

2.用攝動法（perturbation）求近似解析解

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

（為什么要這樣展開？我的理解：把解分成幾部分之和，每一部分的幅值大小由epsilon刻畫，因為epsilon是一個小量，顯然x0(t)對x(t)的貢獻最主要，x1(t)次之，貢獻依次遞減。也就是說認為x(t)主要是按x(t)變化，受到了一些攝動x1(t)，攝動x1(t)對x(t)的影響大小與小量epsilon同階。這種展開應該主要應用在弱非線性的系統，即認為解x(t)仍然主要按對應線性系統的解x0(t)變化，再加上一點非線性的影響（攝動）x1(t),x2(t)，這些攝動的影響是小量的。因為是從小量epsilon對解的影響來看待問題的，所以這種方法主要適用于弱非線性系統）

按這種方式展開之后，為了與原初始條件協調一致，這些分量的初始條件應該為

將展開式（5）式代入原方程（1）式，并讓同一epsilon的冪次的合并，可以得到一系列線性的方程，注意看到由于非線性（這里是三次方）的存在，展開的時候顯然很麻煩，這也是攝動法的一個缺點

這里只寫到epsilon的一次項，一般來說對于近似解這就足夠的，而且更高次的寫出來顯然很麻煩

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

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

對于epsilon的一階項，注意到x0已經求出來了，所以是線性的方程。也就是說用攝動法我們可以得到對每一個epsilon的冪次的線性方程，這是攝動法的優點，把非線性方程分成幾個線性的方程求解。但是要注意到等式右邊一般來說是關于時間t的復雜的表達式（一般是三角函數），所以即使是要求解一個線性微分方程的解析解也是麻煩的

以（7）式子為例，

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

將該結果與上面用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),'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

　　

可見攝動法的解只有在開始的一段很短的時間是對的，很快后面就不對了，會發散（變得很大）。所以攝動法是有問題的。問題出現在什么地方？讓我們看一下用攝動法得到的解

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

secular term之所以會出現是因為在（11）式的右邊（激勵）有這一部分，注意到其激勵頻率等于左邊系統的固有頻率，即發生了共振。

這是攝動法的一個常見缺點（或者說致命缺點？），至于為什么會出現，什么情況下可以避免，這里暫時不深入。總之只需要知道攝動法求解時有時候會由于（人為地引入）共振的出現，導致解有secular term。

3. 平均法

平均法可以用來求解下面這一類關于一階導數非線性的非線性ODE：

上一節用攝動法考慮的例子是一個特例：

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

為什么要這么考慮？因為當非線性消失時，退化的線性ODE的解具有（15）式和（16）式的形式，不過此時的幅值和相位是由初始條件決定的常數；當系統有了一點弱非線性時，我們可以樂觀地假設方程的解仍然具有（15）式的形式，不過此時幅值和相位是會隨時間變化的函數，并且認為幅值和相位是緩慢變化的。為什么認為是緩慢變化的？因為緩慢變化的話方程的解才大致地具有和線性時候的解一樣的形式，只不過由于非線性帶來幅值和相位的緩慢變化。注意這里頻率仍然是固定的。若非線性出現在會影響頻率的剛度項，又該怎么做呢？

這是平均法的核心，或者說亮點吧，就是洞見問題的解具有（15）和（16）式的形式。

將（15）式求導代入（16）式得到一個方程，將（16）式求導代入（1）式，得到兩個耦合的ODEs

整理成

到這一步為止都是精確的，不過是假設了解的形式，但沒有作任何的近似或者截斷。

接下來為了求解（20）和（21）式，要做一點假設：

因為是弱非線性，即epsilon是個小量，所以認為幅值和相位的變化是緩慢的，即幅值和相位的導數是小量。這樣的話，在振動的一個周期內（這里頻率是1，所以周期是2pi），（20）式和（21）式的方程右邊幾乎保持不變是一個常數，這樣就可以用它們在一個周期內的平均作近似。（這樣的近似還是不是很理解）

（20）和（21）式變為下面的形式，是平均化的方程

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

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

最終（23）和（24）式變為

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

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

和數值解作比較發現是較為精確地！而且和數值解相比具有解析形式的優勢。

代碼如下：

% 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. 接下來考慮另一個系統，這是一個著名的非線性系統，叫做van der pol振子，具有（13）式的形式，因此也可以用平均法求近似解析解

求數值解（ode45）時仍然要先轉化為一階的ODEs

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

第一式求導代入第二式，第二式求導代入方程，得到

整理：

接下來做平均：

計算積分同樣要用到三角恒等式

最后得到

求解ODE（實際還很難求解的）：

和數值解作比較，在前面一段時間有點誤差，后面對得還好：

代碼如下

% 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

　　

總結：

平均法的兩個步驟：

一是假設解的形式，這個形式由非線性激勵和/或非線性阻尼中的小參數取零時確定，即系統在無阻尼無激勵時的精確解形式，然后令幅值和相位變為時間的函數

二是認為幅值和相位的變化是小量，在一個周期內不變，所以可以對時間在一個周期上上作平均，去掉快變量

練習：
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.

我的解答

數值解對比

　代碼：

% 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.

看起來有點麻煩，罷了

后續應對平均法作深入了解，是求解非線性微分方程的近似解析解的主要方法

　

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。