[運算] 二階聯立微分方程式

[jimmyjean (色仙)]

我要解一個聯立微分方程的數值解並繪圖
題目是
X''=-0.0675*X'*sqrt(X'^2+Y'^2)
Y''=-9.81-0.0675*Y'*sqrt(X'^2+Y'^2)
初始條件X'(0)=134 X(0)=Y(0)=Y'(0)=0
我降階後的程式如下
function dx=odedx(t,x)
a=0.06725;
g=9.81;
%x(1)=x
%x(2)=x'
%x(3)=y
%x(4)=y'
dx=[x(2);x(4);-a*x(2)*sqrt((x(2))^2+(x(4))^2);-g-a*x(4)*sqrt((x(2))^2+(x(4))^2
)];
end
------------------------------
執行檔如下
x0=[0;0;134;0];
tspan=0:0.025:5;
tol=1e-6;
options=odeset('RelTol',tol,'AbsTol',[tol,tol]);
[t,x]=ode45('odedx',tspan,x0);
plot(x(1),x(3),'-r')
但出來的結果是錯的 所有值幾乎都一樣 但我看不太出來程式碼哪裡有問題 想請教版上
高手 謝謝

--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 140.113.124.28
※ 文章網址: https://www.ptt.cc/bbs/MATLAB/M.1514740513.A.01E.html
推 sunev: dx順序不對
※ 編輯: jimmyjean (140.113.124.28), 01/01/2018 01:53:24
不好意思不是很懂您的說法 是指dx內的順序與初始值順序不一樣嗎
※ 編輯: jimmyjean (140.113.124.28), 01/01/2018 01:54:47
推 sunev: dx(1) = x(2), dx(2) = x(4) ?
我的寫法是 dx=[x';y';x'';y''] 我也有改成dx=[x';x'';y';y'']並修改初始值順序跑過
 但結果是一樣的 不知道您的意思是修改成什麼樣子 謝謝
※ 編輯: jimmyjean (140.113.124.28), 01/01/2018 20:22:35
推 sunev: 如果dx=[x';y';x'';y'']，那x=[x;y;x';y']，dx(1)=x(3)才對
嗯我懂您的意思了 但我改成以下
%x(1)=x
%x(2)=y
%x(3)=x'
%x(4)=y'
dx=[x(3);x(4);-a*x(3)*sqrt((x(3))^2+(x(4))^2);-g-a*x(4)*sqrt((x(3))^2+(x(4))^2
)];
還是無法有一樣的問題 繪出來的圖只有一點 約在(0,3)的位置
※ 編輯: jimmyjean (140.113.124.28), 01/01/2018 21:07:22
推 sunev: plot的指令有問題，你要不要先看看x(1)長什麼樣子
→ LiamIssac: 好習慣 把每個函數都印出來看一下 確認跟手算的都一樣
→ LiamIssac:  再繼續下一步