查看: 6854|回复: 5|关注: 0

[已解决] 关于三元二阶常系数微分方程的求解

[复制链接]

新手

7 麦片

财富积分


050


1

主题

3

帖子

0

最佳答案
求大神帮忙看一下
一个三元二阶常系数微分方程组,我用ode45求得解和某篇论文里面的解析解结果不一致,是哪里出问题了

三元二阶常系数微分方程组

三元二阶常系数微分方程组

通解.png
出处《吴幼明, 罗旗帜. 一类二阶常系数微分方程组的通解[J]. 佛山科学技术学院学报(自然科学版), 2002, 20(2):10-14.
--------------------------------------------------
matlab代码如下


  1. clear
  2. clc

  3. %% 数值解
  4. f=@(t,y)([y(4);
  5.           y(5);
  6.           y(6);
  7.           y(4)-y(5)-y(6)+y(1);
  8.          -y(4)+y(5)-y(6)+y(2);
  9.          -y(4)-y(5)+y(6)+y(3);]);
  10. [t,Y]=ode45(f,[0 2],[1 0 0 0 0 0]);
  11. n=find(Y(:,1)>=0);
  12. plot(t,Y(1:max(n),1),t,Y(1:max(n),2),t,Y(1:max(n),3))
  13. hold on


  14. %% 解析解
  15. a=(1+sqrt(2))/(3*sqrt(2));
  16. b=(sqrt(2)-1)/(3*sqrt(2));
  17. c=(sqrt(5)-1)/(6*sqrt(5));
  18. d=(sqrt(5)+1)/(6*sqrt(5));

  19. a2=(-1-sqrt(2))/(6*sqrt(2));
  20. b2=(1-sqrt(2))/(6*sqrt(2));
  21. c2=(sqrt(5)-1)/(6*sqrt(5));
  22. d2=(sqrt(5)+1)/(6*sqrt(5));

  23. x=linspace(0,2,100);
  24. Y2=a*exp((sqrt(2)-1)*x)+b*exp((-sqrt(2)-1)*x)+c*exp((1+sqrt(5))/2*x)+d*exp((1-sqrt(5))/2*x);
  25. Y3=a2*exp((sqrt(2)-1)*x)+b2*exp((-sqrt(2)-1)*x)+c2*exp((1+sqrt(5))/2*x)+d2*exp((1-sqrt(5))/2*x);
  26. plot(x,Y2,'--',x,Y3,'-.')
  27. legend('x','y','z','x2','y2')

复制代码


结果对比:
结果.png

论坛优秀回答者

专家

2431 麦片

财富积分


20003000


1

主题

2699

帖子

547

最佳答案
  • 关注者: 96
发表于 2018-6-26 19:17:22 | 显示全部楼层
本帖最后由 maple1314168 于 2018-6-26 19:26 编辑

你确认解析解没有错了?
我查你给的论文,好像没有这个例子。(例子的数值解与解析解吻合很好)
是你自己解的?

论坛优秀回答者

中级

799 麦片

财富积分


5001500


2

主题

646

帖子

141

最佳答案
  • 关注者: 46
发表于 2018-6-26 19:20:38 | 显示全部楼层
因为解析解是错的
  1. clear;
  2. f=@(t,y)([y(4);
  3.           y(5);
  4.           y(6);
  5.           y(4)-y(5)-y(6)+y(1);
  6.          -y(4)+y(5)-y(6)+y(2);
  7.          -y(4)-y(5)+y(6)+y(3);]);
  8. [t,Y]=ode45(f,[0 2],[1 0 0 0 0 0]);

  9. x =@(t) (1/30)*(10*exp((1 - sqrt(2))*t) + ...
  10.     5*sqrt(2)*exp((1 - sqrt(2))*t) + 10*exp((1 + sqrt(2))*t) - ...
  11.     5*sqrt(2)*exp((1 + sqrt(2))*t) + ...
  12.     5*exp((1/2)*(-1 - sqrt(5))*t) - ...
  13.     sqrt(5)*exp((1/2)*(-1 - sqrt(5))*t) + ...
  14.     5*exp((1/2)*(-1 + sqrt(5))*t) + ...
  15.     sqrt(5)*exp((1/2)*(-1 + sqrt(5))*t));
  16.      
  17. y =@(t) (-(1/240))*(-1 + sqrt(5))*(1 + sqrt(5))*...
  18.     (10*exp((1 - sqrt(2))*t) + 5*sqrt(2)*exp((1 - sqrt(2))*t) + ...
  19.     10*exp((1 + sqrt(2))*t) - 5*sqrt(2)*exp((1 + sqrt(2))*t) - ...
  20.     10*exp((1/2)*(-1 - sqrt(5))*t) + ...
  21.     2*sqrt(5)*exp((1/2)*(-1 - sqrt(5))*t) - ...
  22.     10*exp((1/2)*(-1 + sqrt(5))*t) - ...
  23.     2*sqrt(5)*exp((1/2)*(-1 + sqrt(5))*t));

  24. z =@(t) (1/60)*(-10*exp((1 - sqrt(2))*t) - ...
  25.     5*sqrt(2)*exp((1 - sqrt(2))*t) - 10*exp((1 + sqrt(2))*t) + ...
  26.     5*sqrt(2)*exp((1 + sqrt(2))*t) + ...
  27.     10*exp((1/2)*(-1 - sqrt(5))*t) - ...
  28.     2*sqrt(5)*exp((1/2)*(-1 - sqrt(5))*t) + ...
  29.     10*exp((1/2)*(-1 + sqrt(5))*t) + ...
  30.     2*sqrt(5)*exp((1/2)*(-1 + sqrt(5))*t));

  31. n=find(Y(:,1)>=0);
  32. subplot(131);plot(t,Y(1:max(n),1),'ks-',t,x(t),'k.-.','MarkerSize',8,'LineWidth',0.5);legend('x','x2')
  33. subplot(132);plot(t,Y(1:max(n),2),'rs-',t,y(t),'r.-.','MarkerSize',8,'LineWidth',0.5);legend('y','y2')
  34. subplot(133);plot(t,Y(1:max(n),3),'bs-',t,z(t),'b.-.','MarkerSize',8,'LineWidth',0.5);legend('z','z2')
复制代码

数值解和正确的解析解两者是重合的

新手

7 麦片

财富积分


050


1

主题

3

帖子

0

最佳答案
 楼主| 发表于 2018-6-26 22:53:03 | 显示全部楼层
maple1314168 发表于 2018-6-26 19:17
你确认解析解没有错了?
我查你给的论文,好像没有这个例子。(例子的数值解与解析解吻合很好)
是你自己解的 ...

解析解的例子是论文里的算例
数值解释我自己写的

新手

7 麦片

财富积分


050


1

主题

3

帖子

0

最佳答案
 楼主| 发表于 2018-6-26 23:14:12 | 显示全部楼层
TouAkira 发表于 2018-6-26 19:20
因为解析解是错的

数值解和正确的解析解两者是重合的

感谢!
请问你这的解析解是怎么得来的
另外,请问你这有没有此类非齐次微分方程组的通解的解析解

论坛优秀回答者

中级

799 麦片

财富积分


5001500


2

主题

646

帖子

141

最佳答案
  • 关注者: 46
发表于 2018-6-26 23:44:49 | 显示全部楼层 |此回复为最佳答案
kercc1985 发表于 2018-6-26 11:14
感谢!
请问你这的解析解是怎么得来的
另外,请问你这有没有此类非齐次微分方程组的通解的解析解 ...

这类不好说,书不在手边没法查,顶楼的ode组是能求的

  1. clear;
  2. syms x(t) y(t) z(t)
  3. D2x = diff(x,t,2);
  4. D2y = diff(y,t,2);
  5. D2z = diff(z,t,2);
  6. Dx = diff(x,t);
  7. Dy = diff(y,t);
  8. Dz = diff(z,t);
  9. eqns = [D2x - Dx + Dy + Dz - x == 0,...
  10.     D2y + Dx - Dy + Dz - y == 0,...
  11.     D2z + Dx + Dy - Dz - z == 0];
  12. cond = [x(0) == 1, y(0) == 0, z(0) == 0, Dx(0) == 0, Dy(0) == 0, Dz(0) == 0];
  13. [x,y,z] = dsolve(eqns,cond);
  14. x=matlabFunction(x);
  15. y=matlabFunction(y);
  16. z=matlabFunction(z);
复制代码

根据版本不同可能需要调整部分函数的调用格式
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

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

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