关于ode45用法的一点浅见和问题

kb24lal

方程：x''+2x'+x=0;满足初始条件(x|t=0)=4、(x'|t=0)=-2于是根据上式可写成以下形式：
function ydot=weifen(t,y)
ydot=zeros(2,1);
ydot(1)=y(2);
ydot(2)=-2*y(2)-y(1); %保存文件


然后就是求解了：
clc;clear
tspan=[0 20];
yo=[-2;0];        %初始条件:这里要用的是(x'|t=0)=-2以及(x''|t=0)=0而不是原条件！此处弄巧成拙，是错误的做法！ode45没这么麻烦，希望大家别被我误导，正确做法请看5#楼！
[t,y]=ode45(@weifen,tspan,yo);
plot(t,y(:,1),t,y(:,2),'--')

————————————————————至此结束
需要说的是，对于和我一样刚开始接触ode45的童鞋，初始条件那块要弄清楚；另外，解出来的y的列向量是x'和x''，而不是x。
这些都是很基础的一些认识，由于我没有搜到讲这些基础的例子，因此自己琢磨了一晚上才搞清楚，浪费了不少时间，因此，以上的这些浅见只适合入门新手 ，大牛可一笑而过～
以上权当抛砖引玉，接下来就是我关心的问题了：
1、我们通常去计算弹簧振子的响应，需要关心位移的时程变化，即相当于上式中的x，那么是不是用ode45就只能得到速度和加速度，位移就得不到了呢？如果可以，请各位大牛指导，谢谢！
2、ode45求得的这些解，有办法让其通过函数的形式表示出来吗，而不是一个个点。

xyz999

对于和我一样刚开始接触ode45的童鞋，初始条件那块要弄清楚；另外，解出来的y的列向量是x'和x''，而不是x。
kb24lal 发表于 2011-3-30 10:41

搞错啦，应该是“另外，解出来的y的列向量是 x 和 x‘ ”

messenger

第1个问题其实2#已经回答了，解出来的是x和x'，其中x就是位移；
第2个问题，因为是数值解，所以不能以函数的形式表示出来。不过，你可以试试用拟合来拟合成某种函数形式。

1、我们通常去计算弹簧振子的响应，需要关心位移的时程变化，即相当于上式中的x，那么是不是用ode45就只能得到速度和加速度，位移就得不到了呢？如果可以，请各位大牛指导，谢谢！
2、ode45求得的这些解，有办法让其通过函数的形式表示出来吗，而不是一个个点。
kb24lal 发表于 2011-3-30 10:41

kb24lal

谢谢你们的关注，3#、4# ，但是对你们的观点不敢苟同我做的这个例子，解析解为：
x=(4+2*t)*exp(-t)
x'=(-2-2*t)*exp(-t)
x''=2*t*exp(-t)
附件是为验证我在1#所作的图，---表示用ode45求出来的x''；红圈表示x'的解析解，红圈里面的蓝线为ode45求解出来的y(:,1),可以看出它们很好的吻合；绿圈为x的解析解，与用ode45求出来的y的列向量不吻合。
其实，你可以自己试一下，这个很简单，我也不敢托大，没有根据我当然不会在这里误导他人。

xyz999

方程（和楼主的一样）：
function ydot=weifen(t,y)
ydot=zeros(2,1);
ydot(1)=y(2);
ydot(2)=-2*y(2)-y(1);

求解：
clc;clear
tspan=[0 20];
yo=[4;-2];   % % 初始条件(x|t=0)=4、(x'|t=0)=-2
[t,y]=ode45(@weifen,tspan,yo);
xa = (4+2*t).*exp(-t); % 解析解
plot(t,y(:,1),t,y(:,2),'--',t,xa,'o')
legend('numerical x', 'numerical dx/dt', 'analytical solution x')

结果：

kb24lal

5# xyz999 真的是非常感谢！！！
看来我是对那个初始条件混淆了，在这里献丑了。。。
请问，这个yo给的不是function里面的ydot吗？ydot(1)不就是x'吗？ydot(2)不就是x''吗？

xyz999

yo 不是function ydot=weifen(t,y) 中的ydot，是 0 时刻的 [位移；速度]。ydot(1)是x‘，ydot(2)是x".

kb24lal

yo 不是function ydot=weifen(t,y) 中的ydot，是 0 时刻的 [位移；速度]。ydot(1)是x‘，ydot(2)是x".
xyz999 发表于 2011-3-31 18:15

原来如此，谢谢了！
是不是对于2阶的微分方程，ode45所需要的初始条件yo指的就是x以及x' ；而对于1阶的则是x ？

nwcwww

本科读大2的时候有非线性动力学和混沌这门课，当时基本全靠ode45打天下。。。
MATLAB解ODE的方法相当多，数值和解析解都不在话下。当然这和你的toolbox license也有关。
ode45因为是数值解法，所以无法直接让结果通过函数形式表达出来，楼主可以换个方法直接求解析解。

解析解方面有两种方法我用的多些，optimization toolbox里的fsolve，以及symbolic maths里进mupad直接算。
fsolve解ODE如下：

1. >> sol = dsolve('D2x+2*Dx+x(t)=0','x(0) = 4','Dx(0) = -2','t')
   
2. 
   
3. sol =
   
4. 
   
5. 4/exp(t) + (2*t)/exp(t)

复制代码

Mupad直接截图好了，代码是其中红色的部分（看不清的话可以单击图片放大）：

xyz999

8# kb24lal
对！
其实，对ode45来说，没有“2阶”的概念，都是一阶的微分方程或方程组。为了利用ode45，人们把二阶或二阶以上的方程化成了等效的一阶方程组。

kb24lal

10# xyz999 嗯，谢谢！
由于对2阶的这种处理方法，所以让人感觉ode45很神奇，好像能自动识别似的。。多亏你的耐心解释，让我获益匪浅！

kb24lal

9# nwcwww 嗯，对于有解析解的微分方程这个的确可以，对于大部分非线性微分方程还是只能用数值解法。。3q all the same

南边的北边

用数值方法求积ODE方程时，如果采用ode45的话，必定要将其写为状态空间，如果原系统有n个自由度，那么状态空间是与2*n个方程，当求积的方程数目不多的收，计算量应该不是问题，但如果自由度数很大的，再写成状态方程的话就会增加很大的工作了。
对于一般的大型的微分方程组，还是采用newmark等一类隐式解法吧，

bobohan

如果要把ydot(2),也即X''的数值解画出来，那怎么实现？因为ODE45好像只能求解X, X'
多谢了！ 7# xyz999

kb24lal

bobohan 发表于 2011-6-28 18:14
如果要把ydot(2),也即X''的数值解画出来，那怎么实现？因为ODE45好像只能求解X, X'
多谢了！ 7# xyz999  ...

呵呵，那就直接看我一楼做的结果呀～
讲明白点吧：
x''+2x'+x=0   满足初始条件(x|t=0)=4、(x'|t=0)=-2
令y=x',并对方程两边求导得到：
y''+2y'+y=0   满足初始条件(y|t=0)=-2、(y'|t=0)=0
使用ode45可求出y以及y'，即原方程的x'和x''