微分方程的框架

導彈引例問題

基本概念和怎樣建立微分方程

MATLAB求微分方程解析解

%% 例1
clear;clc
dsolve('y-Dy=2*x','x')  % 這裡要指定自變數為x
% 2*x + C1*exp(x) + 2  (這裡的C1表示任意常數，有時候也會出現C2 C3等)
dsolve('y-Dy=2*x')  % 如果不指定自變數的話，會預設自變數為t，x會看成一個常數
% 2*x + C2*exp(t)

% 注意：最新版本的matlab會逐漸淘汰上面那種寫法（雖然我個人覺得上面的寫法更方便）
% 下面這種寫法是新版的matlab推薦的方式（和我們上一講符號運算中解方程的寫法類似）
syms y(x)
eqn = (y - diff(y,x) == 2*x);    % 注意原來方程中的“=”一定要改成“==”
dsolve(eqn)

%% 如果微分方程中還有其他的未知引數怎麼辦？
% 方法1
dsolve('y-Dy=a*x','x')  % a是一個未知的引數
% 方法2
syms y(x) a
eqn = (y - diff(y,x) == a*x);  
dsolve(eqn)

%% 例2 
% 方法1
dsolve('y-Dy=2*x','y(0)=3','x')
% 2*x + exp(x) + 2
% 方法2
syms y(x)
eqn = (y - diff(y,x) == 2*x);  
cond = (y(0) == 3);
dsolve(eqn,cond)
% 2*x + exp(x) + 2


%% 例3
% 方法1
dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x') 
% 3*sin(5*x)*exp(-2*x)
% 方法2
syms y(x)
eqn = (diff(y,x,2) + 4 *diff(y,x) + 29*y  == 0);  
Dy = diff(y,x); % 定義變數Dy為y的一階導數
cond = [(y(0) == 0) ,(Dy(0) ==15)] ; % 有兩個條件，可以寫到一個向量中儲存
dsolve(eqn,cond)
% 3*sin(5*x)*exp(-2*x)

%% 例4
% 方法1
[x,y,z] = dsolve('Dx=2*x-3*y+3*z+t','Dy=4*x-5*y+3*z+t','Dz=4*x-4*y+2*z+t','t') 

% 方法2
syms x(t) y(t) z(t)
eqn1 = (diff(x,t)  == 2*x-3*y+3*z+t); 
eqn2 = (diff(y,t)  == 4*x-5*y+3*z+t); 
eqn3 = (diff(z,t)  == 4*x-4*y+2*z+t); 
eqns = [eqn1 eqn2 eqn3];
[x,y,z] = dsolve(eqns)
% x = exp(2*t)*(C2- (exp(-2*t)*(2*t + 1))/4) + C3*exp(-t)
% y = exp(2*t)*(C2 - (exp(-2*t)*(2*t + 1))/4) + C3*exp(-t) + C4*exp(-2*t)
% z = exp(2*t)*(C2 - (exp(-2*t)*(2*t + 1))/4) + C4*exp(-2*t)
mupad  % 最新版本matlab可能會報錯，將計算結果複製到裡面，使結果可讀。
% 如果新版matlab用不了mupad的話，可以使用更新13中介紹到的實時指令碼
simplify(y)  % simplify函式可以簡化表示式
latex(y) % 轉換成latex程式碼，複製到Axmath或者word自帶的公式編輯器（低版本不知道支不支援）
% 如果太過於複雜的話可能會報錯，大家可以自己測試

%% 不是所有的微分方程都可以，導彈追擊那一題就沒有解析解
% 假設 v=100
[x,y] = dsolve('Dx = 3*100*(20+sqrt(2)/2*100*t-x)/sqrt((20+sqrt(2)/2*100*t-x)^2+(sqrt(2)/2*100*t-y)^2)','Dy = 3*100*(sqrt(2)/2*100*t-y)/sqrt((20+sqrt(2)/2*100*t-x)^2+(sqrt(2)/2*100*t-y)^2)','x(0)=0,y(0)=0','t')  
% 警告: Explicit solution could not be found.