微分方程的基本概念
微分方程：一般的，凡表示未知函數、未知函數的導數與自變量之間的關系的方程，叫做微分方程，有時也簡稱方程。
微分方程的階：微分方程中所出現的未知函數的最高階導數的階數，叫做微分方程的階。
微分方程的解：使得微分方程成立的函數稱為微分方程的解。不含任意常數的解稱為微分方程的特解。若微分方程的解中所含的互相獨立的任意常數的個數與微分方程的階數相等，稱這個解為微分方程的通解。
微分方程的通解：如果微分方程的解中含有任意常數，且任意常數的個數與微分方程的階數相同，這樣的解叫做微分方程的通解。
微分方程的解析解
一般隻有比較簡單的微分方程才能求出解析解，後面做模組化一般都是求數值解。
這裡我們講的基本都是直接用MATLAB程式算的
MATLAB求解指令為
dsolve(‘方程1’, ‘方程2’,…‘方程n’, ‘初始條件’, ‘自變量’)
記号:在表達微分方程時，用字母D表示求微分，D2、D3等表示求高階微分.任何D後所跟的字母為因變量，自變量可以指定或由系統規則標明為确省.
例一：求
dudt=1+u2\frac{du}{dt}=1+u^2
dtdu=1+u2通解
dsolve('Du=1+u^2','t')
tan(C1 + t)%通解
1i%虛數+i、-i數值解
-1i
例二：
{d2ydx2+4dydx+29y=0y(0)=0y′(0)=15\begin{cases}
\frac{d^2y}{dx^2}+4\frac{dy}{dx}+29y=0
\\
y(0)=0\\
y'(0)=15
\end{cases}
⎩⎪⎨⎪⎧dx2d2y+4dxdy+29y=0y(0)=0y′(0)=15
y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x')
y =
3*sin(5*x)*exp(-2*x)
[x,y,z]=dsolve('Dx=2*x-3*y+3*z','Dy=4*x-5*y+3*z','Dz=4*x-4*y+2*z', 't')；
x=simplify(x) % 将x化簡
y=simplify(y)
z=simplify(z)
x =
C1*exp(-t) + C3*exp(2*t)
y =
exp(-2*t)*(C2 + C1*exp(t) + C3*exp(4*t))
z =
C2*exp(-2*t) + C3*exp(2*t)
微分方程的數值解
定義
在生産和科研中所處理的微分方程往往很複雜且大多得不出一般解。而在實際上對初值問題，一般是要求得到解在若幹個點上滿足規定精确度的近似值，或者得到一個滿足精确度要求的便于計算的表達式。是以，研究常微分方程的數值解法是十分必要的。
建立數值解法的一些途徑
1、歐拉法。
用差商代替導數
若步長h較小，則有
故有公式：
改進的歐拉法
使用數值積分
對方程 y’=f(x, y), 兩邊由xi到xi+1積分，并利用梯形公式，有：
實際應用時，與歐拉公式結合使用：
使用泰勒公式
以此方法為基礎，有龍格-庫塔法、線性多步法等方法。具體的太麻煩了，用的時候直接用MATLAB調函數就行。k越大，則數值公式的精度越高。
數值公式的精度
當一個數值公式的截斷誤差可表示為$O(h^{k+1})時(k為正整數，h為步長)，稱它是一個k階公式。
歐拉法是一階公式，改進的歐拉法是二階公式。
龍格-庫塔法有二階公式和四階公式。
線性多步法有四階阿達姆斯外插公式和内插公式。
用Matlab軟體求常微分方程的數值解
具體詳細的可以help一下看看，我也不是多明白，嘿嘿，就不說了。
注意：
在解n個未知函數的方程組時，x0和x均為n維向量，m-檔案中的待解方程組應以x的分量形式寫成.
使用Matlab軟體求數值解時，高階微分方程必須等價地變換成一階微分方程組.
沒看懂？不要緊，看個例子就會了
{d2xdt2−1000(1−x2)dxdt+x=0x(0)=2x′(0)=0\begin{cases}
\frac{d^2x}{dt^2}-1000(1-x^2)\frac{dx}{dt}+x=0
\\
x(0)=2\\
x'(0)=0
\end{cases}
⎩⎪⎨⎪⎧dt2d2x−1000(1−x2)dtdx+x=0x(0)=2x′(0)=0
由上面第一個式子可以知道，這是一個二階微分方程，我們先對其降階。
令
y1=x,y2=x′y_1=x,y_2=x'
y1=x,y2=x′ 然後帶入上面的方程可以得到
MATLAB求解過程：
建立 .m 檔案vdp1000.m如下：
function dy=vdp1000(t,y)
dy=zeros(2,1);%初始化
dy(1)=y(2);
dy(2)=1000*(1-y(1)^2)*y(2)-y(1);
2、取t0=0，tf=3000，輸入指令：
[T,Y]=ode15s('vdp1000',[0 3000],[2 0]); %[2 0]表示初值
plot(T,Y(:,1),'-')
結果如圖所示，具體T Y的值看工作區
模組化執行個體解析
飛彈追蹤問題
設位于坐标原點的甲艦向位于x軸上點A(1, 0)處的乙艦發射飛彈，飛彈頭始終對準乙艦.如果乙艦以最大的速度v0(是常數)沿平行于y軸的直線行駛，飛彈的速度是5v0，求飛彈運作的曲線方程.又乙艦行駛多遠時，飛彈将它擊中？
解法一、解析法
設t時刻，飛彈的坐标為
P:(x(t),y(t)),P:(x(t),y(t)),
P:(x(t),y(t)),乙艦坐标為
Q:(1,v0t)Q:(1,v_0t)
Q:(1,v0t)示意圖如下
可以得到
y′y'
y′的導數即
tan(θ)tan(\theta)
tan(θ)為：
y′=tan(θ)=v0t−y1−xy'=tan(\theta)=\frac{v_0t-y}{1-x}
y′=tan(θ)=1−xv0t−y (1)
化簡為
v0t=(1−x)y′+yv_0t=(1-x)y'+y
v0t=(1−x)y′+y
根據題意，飛彈的速度是戰艦的5倍，是以飛彈的弧長，為戰艦路程的5倍(時間一樣都為t)，由弧長公式可以得到：
聯立公式1，2，對x求導(注意這裡不是求偏導，y是x的函數，是以
((1−x)y′)′=−y′+(1−x)y′′((1-x)y')'=-y'+(1-x)y''
((1−x)y′)′=−y′+(1−x)y′′，)可以化簡為
初始條件為
y(0)=0,y′(0)=0y(0)=0,y'(0)=0
y(0)=0,y′(0)=0
這裡就可以求y的解析式了，也就是這個微分方程的通解
MATLAB代碼和結果為
y=dsolve('(1-x)*D2y-sqrt(1+Dy^2)/5=0','y(0)=0','Dy(0)=0','x')
y=simplify(y)
(5*(-1)^(1/5)*(x - 1)^(4/5))/8 + (5*(-1)^(4/5)*(x - 1)^(6/5))/12 + 5/24
- (5*(-1)^(1/5)*(x - 1)^(4/5))/8 - (5*(-1)^(4/5)*(x - 1)^(6/5))/12 - 5/24
可以看到，有兩個結果，上面那個是沿着y軸正向行駛，下面是負向行駛
當兩者的坐标相等時，兩者相遇，是以x=1時
syms x
f(x)= (5*(-1)^(1/5)*(x - 1)^(4/5))/8 + (5*(-1)^(4/5)*(x - 1)^(6/5))/12 + 5/24;
f(1)
ans =
5/24
是以擊中時的坐标為(1，5/24)，也就是行駛5/24時被擊中
解法二、數值解
令y1=y,y2=y1’，将方程(3)化為一階微分方程組。
然後MATLAB求數值解，先建立.m函數
function dy=eq1(x,y)
dy=zeros(2,1);
dy(1)=y(2);
dy(2)=1/5*sqrt(1+y(1)^2)/(1-x);
再求解畫圖
[x,y]=ode15s('eq1',[0 1],[0 0]);
plot(x,y(:,1),'b')
hold on
y=0:0.1:2;
plot(1,y,'b^')
y
2.0000
最後一個y值為0.2000與5/24=0.2083相差不大
解法三、建立參數方程求數值解
設時刻t乙艦的坐标為(X(t)，Y(t))，飛彈的坐标為(x(t)，y(t)).設飛彈速度恒為w,則 (dxdt)2+(dydt)2=w2(\frac{dx}{dt})^2+(\frac{dy}{dt})^2=w^2(dtdx)2+(dtdy)2=w2
由于彈頭始終對準乙艦，故飛彈的速度平行于乙艦與飛彈頭位置的差向量，二者成比例
将(2)分開算帶入(1)，消去參數得
因乙艦以速度v0沿直線x=1運動，設v0=1，則w=5，X=1，Y=t
是以飛彈運動軌迹的參數方程為：
令x=yx,y=y2x=y_x,y=y_2x=yx,y=y2
%.m函數檔案來
function dy=eq2(t,y)
dy=zeros(2,1);
dy(1)=5*(1-y(1))/sqrt((1-y(1))^2+(t-y(2))^2);
dy(2)=5*(t-y(2))/sqrt((1-y(1))^2+(t-y(2))^2)
%主程式
[t,y]=ode45('eq2',[0 2],[0 0]);
Y=0:0.01:2;
plot(1,Y,'-'), hold on
plot(y(:,1),y(:,2),'*')
1.9891
結果 1.9891 三種算法結果都差不多