作者jimmyjean (色仙)
看板MATLAB
标题[运算] 二阶联立微分方程式
时间Mon Jan 1 01:15:11 2018
我要解一个联立微分方程的数值解并绘图
题目是
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://webptt.com/cn.aspx?n=bbs/MATLAB/M.1514740513.A.01E.html
1F：推 sunev: dx顺序不对01/01 01:38
※ 编辑: jimmyjean (140.113.124.28), 01/01/2018 01:53:24
不好意思不是很懂您的说法 是指dx内的顺序与初始值顺序不一样吗
※ 编辑: jimmyjean (140.113.124.28), 01/01/2018 01:54:47
2F：推 sunev: dx(1) = x(2), dx(2) = x(4) ?01/01 15:59
我的写法是 dx=[x';y';x'';y''] 我也有改成dx=[x';x'';y';y'']并修改初始值顺序跑过
但结果是一样的 不知道您的意思是修改成什麽样子 谢谢
※ 编辑: jimmyjean (140.113.124.28), 01/01/2018 20:22:35
3F：推 sunev: 如果dx=[x';y';x'';y'']，那x=[x;y;x';y']，dx(1)=x(3)才对01/01 20:35
嗯我懂您的意思了 但我改成以下
%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
4F：推 sunev: plot的指令有问题，你要不要先看看x(1)长什麽样子 01/01 23:28
5F：→ LiamIssac: 好习惯 把每个函数都印出来看一下 确认跟手算的都一样 01/01 23:40
6F：→ LiamIssac: 再继续下一步 01/01 23:40