查看: 124|回复: 8|关注: 0

[已解决] 用dsolve求解微分方程组遇到问题

[复制链接]

新手

10 麦片

财富积分


050


5

主题

17

帖子

0

最佳答案
qm = 1; qe = 0.2; b0 = 1; b2 = -1; k0 = 1;k2 = 0; xm = 0.8; r = 1/qm+1/qe; r2 = 2*b2/qe; r4 = b2^2/qe; qt = 1/r;
k = @(x)(x+k2*x.^3);   b = @(x)(1+b2*x.^2);   r = @(x)(1/qm+(1+b2*x.^2).^2/qe);

y0 = @(z)1/(4-8*z^2+4*r0^2*z^2+4*z^4);
z = 1; T = 2*pi/z; Ts = T/48; Tmax = 5*T; zm2 = 1/qt^2+(z-1/z)^2; x00 = -(1-1/z^2)/zm2; v0 = 1/qt/zm2;

%%%%%%%%%%%%%%以上都是一些变量,函数句柄的定义
问题是,想求解下述微分方程组,自变量是时间t,上面 k , b , r 的自变量x是位移,x也是时间t的函数

dx/dt=v;                                                              %位移对时间求导,是速度
dv/dt=-r(x(t))*v(t)-k(x(t))+b(x(t))*cos(z*t);       %速度对时间求导,是加速度
初始条件:x(0)=x00, x(0)=v0      
{x,v},{t,0,Tmax}

如何用dsolve求解x, v, a,即要得到位移,速度和加速度,在时间t 从0到Tmax的值
求大神指导!

论坛优秀回答者

5

主题

1918

帖子

552

最佳答案
  • 关注者: 159
发表于 2020-3-23 10:55:06 | 显示全部楼层
回去检查吧:
r = 1/qm+1/qe;
r = @(x)(1/qm+(1+b2*x.^2).^2/qe);
这肯定不行啊,变量跟函数重名会直接导致不可预测的结果的。
提问请:①准确描述问题②提出你的思考(等着抄作业的一律锁帖)③提供代码文本而非截图④及时反馈

新手

10 麦片

财富积分


050


5

主题

17

帖子

0

最佳答案
 楼主| 发表于 2020-3-23 10:59:12 | 显示全部楼层
TouAkira 发表于 2020-3-23 10:55
回去检查吧:
r = 1/qm+1/qe;
r = @(x)(1/qm+(1+b2*x.^2).^2/qe);

原本是这样的:kx = @(x)(x+k2*x.^3);   bx = @(x)(1+b2*x.^2);   rx = @(x)(1/qm+(1+b2*x.^2).^2/qe);
因为我怕下面dv/dt的表达式看着乱,所以贴出来的代码把kx rx 里的 x删掉了

论坛优秀回答者

5

主题

1918

帖子

552

最佳答案
  • 关注者: 159
发表于 2020-3-23 11:22:43 | 显示全部楼层 |此回复为最佳答案
guowy 发表于 2020-3-22 22:59
原本是这样的:kx = @(x)(x+k2*x.^3);   bx = @(x)(1+b2*x.^2);   rx = @(x)(1/qm+(1+b2*x.^2).^2/qe);
...

1.不要随便删东西,特别是在删掉后会混淆的时候
2.dsolve一般只用来求解析解的,这种复杂的微分方程组一般要用ode solver求数值解,常用的是ode45
www.mathworks.com/help/matlab/ref/ode45.html
大致就是下面这样的(示意,具体细节自己照着文档的示例改)
  1. Myfun = @( t, y ) [ y( 2 );
  2.     -rx( y( 1 ) ) * y( 2 ) - kx( y( 1 ) ) + bx( y( 1 ) ) * cos( z * t ) ];
  3. [ T, sol ] = ode45( Myfun, [ 0, Tmax ], [ x00, v0 ] );
  4. plot( T, sol( :, 1 ), 'k-', T, sol( :, 2 ), 'r-' )
复制代码
提问请:①准确描述问题②提出你的思考(等着抄作业的一律锁帖)③提供代码文本而非截图④及时反馈

新手

10 麦片

财富积分


050


5

主题

17

帖子

0

最佳答案
 楼主| 发表于 2020-3-23 11:48:26 | 显示全部楼层
TouAkira 发表于 2020-3-23 11:22
1.不要随便删东西,特别是在删掉后会混淆的时候
2.dsolve一般只用来求解析解的,这种复杂的微分方程组一 ...

感谢!
主要是 Myfun = @( t, y ) [ y( 2 );-rx( y( 1 ) ) * y( 2 ) - kx( y( 1 ) ) + bx( y( 1 ) ) * cos( z * t ) ]; 这句不会写
尝试写过,但总是与您的答案这里差一些,那里有点不同
看了您的答案有种“原来不是很复杂的”的感觉,是我想复杂了

新手

10 麦片

财富积分


050


5

主题

17

帖子

0

最佳答案
 楼主| 发表于 2020-3-23 15:40:02 | 显示全部楼层
TouAkira 发表于 2020-3-23 11:22
1.不要随便删东西,特别是在删掉后会混淆的时候
2.dsolve一般只用来求解析解的,这种复杂的微分方程组一 ...

版主,我又来了:lol
上面通过求微分方程组求得x(t)的数值解,也就是sol(:,1),现在我要对x(t)进行数值积分,在(Tmax-T,Tmax)范围内
integral()的参数需要是函数句柄;cumtrapz()吧,mathworks里的例子显示结果与真值有差异,

                               
登录/注册后可看大图

还是说对sol(:,1)进行拟合,定义函数句柄,再用integral()?

新手

10 麦片

财富积分


050


5

主题

17

帖子

0

最佳答案
 楼主| 发表于 2020-3-23 16:06:25 | 显示全部楼层
TouAkira 发表于 2020-3-23 11:22
1.不要随便删东西,特别是在删掉后会混淆的时候
2.dsolve一般只用来求解析解的,这种复杂的微分方程组一 ...

版主,我又来了~~
上面对x(t)通过微分方程求得数值解,即sol(:,1),现在要对x(t)进行数值积分,范围(Tmax-T, T)
integral()第一个参数要是函数句柄,cumtrapz()得到的结果和真值有一定差异,这个要怎么解决?

新手

10 麦片

财富积分


050


5

主题

17

帖子

0

最佳答案
 楼主| 发表于 2020-3-23 16:40:15 | 显示全部楼层
TouAkira 发表于 2020-3-23 11:22
1.不要随便删东西,特别是在删掉后会混淆的时候
2.dsolve一般只用来求解析解的,这种复杂的微分方程组一 ...

用cumtrapz()还有一个问题,就是我积分范围是(Tmax-T, Tmax) ,上面例子的结果sol(:,1)里面没有对应的 Tmax-T 下的值
求指导!

新手

10 麦片

财富积分


050


5

主题

17

帖子

0

最佳答案
 楼主| 发表于 2020-3-23 17:20:30 | 显示全部楼层
TouAkira 发表于 2020-3-23 11:22
1.不要随便删东西,特别是在删掉后会混淆的时候
2.dsolve一般只用来求解析解的,这种复杂的微分方程组一 ...

能不能把求得的sol(:,1)转为函数句柄,后面需要
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

站长推荐上一条 /4 下一条

快速回复 返回顶部 返回列表