之前過冷水有和大家分享熱傳導方程求解的方法，其本質上是微分方程的問題。考慮大多數讀者對微分方程求解方法比較陌生，所以過冷水本期簡單普及一下微分方程的求解問題。
關於微分方程你需要了解：含有未知的函數及其某些階的導數以及其自變量本身的方程稱為微分方程。如果未知函數是一元函數，則稱為常微分方程。如果未知函數是多元函數，則稱為偏微分方程。聯繫一些未知函數的一組微分方程稱為微分方程組。微分方程中出現的未知函數的導數的最高階稱為微分方程的階。
有些微分方程比較簡單可直接通過積分求解。例如一階常係數線性常微分方程：
syms y(x) a beqn = diff(y,x) == a*y+b;S = dsolve(eqn)S =-(b - C3*exp(a*x))/asyms y(x) a beqn = diff(y,x,2)+2*diff(y,x,1)+4*y ==0；S = dsolve(eqn)eqn(x) =4*y(x) + 2*diff(y(x), x) + diff(y(x), x, x) == 0S =C4*exp(-x)*cos(3^(1/2)*x) + C5*exp(-x)*sin(3^(1/2)*x)演示了兩個比較簡單的微分方程用符號解微分方程的方法解出通解，在我們實際問題中少數特殊方程可用初等積分法求解外，大部分微分方程無顯示解，應用主要依靠數值解法。考慮一階常微分方程組初值問題：
其中y=(y1,y2,...,ym)T，f=(f1,f2,...,fm)T，y0=(y10,y20,...,ym0)T，所謂數值解，就是尋求解y(t)在一些列離散節點t0<t1<t2<...<tn<tf 上的近似值yk(k=0,1,...,n).稱hk=tk+1-tk為步長，已知：
求其數值解。自己根據差分方程思想編程如下：
clear allwarning offfeature jit offf=inline('y-2*x/y','x','y');a=0;b=1;h=0.1;n=(b-a)./h;x=zeros(1,n+1);y=zeros(1,n+1);y(1)=1;for i=1:n+1 x(i)=a+(i-1)*h; y(i+1)=y(i)+h*f(x(i),y(i));endy=y(1:n+1);一般來講符號法的運算會比單純的數值運算可具有科學準確性。因為該問題比較簡單，可以採用符號微分法求解，用符號計算為對比看差分法數值運算精度如何。代碼如下：
%符號法解微分方程y1=dsolve('Dy=y-2*x/y','y(0)=1','x')%繪圖對比figure1 = figure;axes1 = axes('Parent',figure1);hold(axes1,'on');plot(x,y,'DisplayName','差分法','LineWidth',2);h1=ezplot(y1,[0,1]);set(h1,'DisplayName','符號法','LineWidth',2);xlabel('$x$','FontWeight','bold','Interpreter','latex');ylabel('$f(x)$','Interpreter','latex');title('差分法&符號法 解微分方程對比');xlim(axes1,[0 1]);ylim(axes1,[0.93232679283255 1.79972401473633]);set(axes1,'FontSize',14,'FontWeight','bold');legend1 = legend(axes1,'show');set(legend1,'Position',[0.148120235966763 0.82820510022613 0.127868850472195 0.0933333308498066]);我們再來看一個案例：
擬採用兩種符號運算方法，兩種數值運算方法。代碼如下：%第一種解微分方程的方法syms y(x) a beqn1 = diff(y,x,1)==-2*y(x)+2*x^2+2*x;eqn2 = y(0)==1;y1= dsolve([eqn1,eqn2]);%第二種解微分方程的方法y2=dsolve('Dy=-2*y+2*x^2+2*x','y(0)=1','x');%第三種解微分方程的方法warning offfeature jit offf=inline('-2*y+2*x^2+2*x','x','y');a=0;b=0.5;h=0.0001;n=(b-a)./h;x=zeros(1,n+1);y=zeros(1,n+1);y(1)=1;for i=1:n+1; x3(i)=a+(i-1)*h; y(i+1)=y(i)+h*f(x(i),y(i));endy3=y(1:n+1);%第四種微分方程數值解[x4,y4]=ode23(f,[0 0.5],1);%繪圖figure1 = figure;axes1 = axes('Parent',figure1);hold(axes1,'on');h1=ezplot(y1,[0,0.5]);h2=ezplot(y2,[0,0.5]);set(h1,'DisplayName','$y_1 =exp(-2x) + x^2$','LineWidth',2);set(h2,'DisplayName','$y_2=2x + 2exp(-2x) - x^2 - 1$','LineWidth',2);plot(x3,y3,'DisplayName','$(x_3,y_3)$','LineWidth',2);plot(x4,y4,'DisplayName','$y_4=ode23(f(x))$','MarkerFaceColor',[0.494117647409439 0.184313729405403 0.556862771511078],'MarkerSize',8,'Marker','o','LineStyle','none');xlabel('$x$','FontSize',16,'Interpreter','latex');ylabel('$y_4$','FontSize',20,'Interpreter','latex');title('微分方程不同解法數值解對比圖');xlim(axes1,[0 0.5]);ylim(axes1,[0.43978341778928 1.0459754645536]);set(axes1,'FontSize',14,'LineWidth',2);legend1 = legend(axes1,'show');set(legend1,'Position',[0.604420982252212 0.73205553519439 0.220018931648527 0.166277798138944],'Interpreter','latex','FontSize',14);本期推文過冷水就是想講一下簡單的微分方程求解的方法，讓大家足夠解決常見問題就OK了！至於那些複雜問題，萬丈高樓平地起，Monte Carlo算法不也講了好幾期的嗎？敬請期待下期的複雜偏微分方程組的求解方法。
往期回顧>>>>>>歡迎各路英雄豪傑來搞
積分變量替換到legendre微分變換
數值計算——MATLAB數值積分原理詳講
數值優化—三種複雜函數數值積分方法實例演示
非線性方程組求解迭代算法&圖像尋初始值講解
在matlab愛好者公眾號中，回復「QQ」加入公眾號專屬Q群；回復「原創」獲取小編原創代碼；回復「星球」加入資源分享園地知識星球。
如需轉載，請在公眾號中回復「轉載」獲取授權，未經授權擅自搬運抄襲的，必將追究其責任！