声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2289|回复: 11

[编程技巧] 高手救急@求解变参数常微分方程

[复制链接]
发表于 2009-3-26 19:54 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
问题如下,在MATLAB环境里解常微分方程组

dx/dt=a*x+b
dy/dt=c*x
求t∈[0,10]的数值解,步长为1
a,b,c为参数,且方程运算的每一步参数的值均有变化。参数a的变化有规律,b,c的变化没有规律,必须从外部输入。参数和初值如下:
t       a       b     c       x     y
0     1.0     2     1      1     2
1     1.1     6     4
2     1.2     3     2
3     1.3     5     2
4     1.4     3     1
5     1.5     7     3
6     1.6     2     1
7     1.7     3     2
8     1.8     5     4
9     1.9     4     2
10   2.0     2     3
回复
分享到:

使用道具 举报

 楼主| 发表于 2009-3-26 21:12 | 显示全部楼层

补充

这种方程用ODE45能解吗,如果不能的话又用什么方法呢?
期待高手指教
发表于 2009-3-26 21:53 | 显示全部楼层
建議樓主看下本版規則!
 楼主| 发表于 2009-3-26 23:24 | 显示全部楼层
楼上说的在理,正在学习论坛版规,虚心向各位请教:handshake
发表于 2009-3-27 14:35 | 显示全部楼层

  1. function odesss
  2. clear;clc
  3. [T,Y] = ode45(@xdot222,[0 10],[1 2]);
  4. plot(T,Y)
  5. function xd=xdot222(tt,x)
  6. xd=zeros(2,1);
  7. t=0:10;
  8. a=1.0:0.1:2.0;
  9. b =[ 2     6     3     5     3     7     2     3     5     4     2];
  10. c =[1     4     2     2     1     3     1     2     4     2     3];
  11. at=interp1(t,a,tt);
  12. bt=interp1(t,b,tt);
  13. ct=interp1(t,c,tt);
  14. xd(1)=at*x(1)+bt;
  15. xd(2)=ct*x(1);
复制代码
untitled.jpg

评分

1

查看全部评分

 楼主| 发表于 2009-3-27 17:27 | 显示全部楼层

谢谢

谢谢四楼sogooda朋友的热心帮助。正在学习你的方法,我课题里面出现的变参数有十多个,有问题再向你请教,谢谢你,你让我看到了解决这个问题的希望。
 楼主| 发表于 2009-3-27 21:39 | 显示全部楼层

困惑

刚才我计算了一下这个方程的解析解,发现跟四楼的方法得到的数值解相差很大,计算数值解时我把xdot222放到m文件里面,
输入的命令是[T,Y] = ode45(@xdot222,[0:10],[1 2])
两个计算结果附在下面了
自变量     参数      解析解   数值解
tabc    x    yxy
0121    1     200
11.16413.935980727.2217481010
21.23236.0811173350.468528894070
31.352235.5657149339.7934075140230
41.431847.7687091598.2633636560660
51.57310240.9070220411.8140323903460
61.62133219.5085220756.067831129013440
71.732407147.1404478972.75345881060110
81.8546777613.0315061273.18338500578300
91.94282896088.8187259003.921526603211940
10223970330389.814554955551513061020734240

期待高手解答
发表于 2009-3-27 21:59 | 显示全部楼层

回复 7楼 ljunw003 的帖子

先看看趋势对不对,如果只是精度不好,可以通过odeset来改善。
doc odeset
 楼主| 发表于 2009-3-28 09:45 | 显示全部楼层

怎样设置options

趋势是对的,跟解析解有相似的趋势
应该是自己水平太菜的缘故,按照sogooda朋友提供了思路我修改了半天精度还是没有大的改善,期待高手指点一下options怎样设置才使数值解跟解析解最接近。
 楼主| 发表于 2009-3-30 19:24 | 显示全部楼层

问题解决了,回来汇报

感谢大家对我的提问的关注,特别要感谢sogooda 朋友的热心帮助,现在我已经有这个问题的解决办法了,特意回来发到这里,跟大家交流

function dy=bcs(t,y)
dy=zeros(2,1)
global k
dy(1)=k(1)*y(1)+k(2);
dy(2)=k(3)*y(1);

clear
yy=[];y0=[1 2];
global k
for i=1:10
canshu=importdata('can.txt');k=canshu(i,:);
[t,y]=ode45(@bcs,[0:1:1],y0);
yy=[yy;y(end,:)];
y0=y(end,:);
end
save result.txt yy -ascii

不知道什么原因。图片传不上去,老提示格式不对,可是我明明用的是JPG,遗憾

can.txt

95 Bytes, 下载次数: 13

参数

result.txt

340 Bytes, 下载次数: 11

结果

评分

1

查看全部评分

发表于 2011-1-8 10:25 | 显示全部楼层
高手,我也遇到这样一个问题,是关于曲轴连杆动力学方程的,因为压力p也是个变参数,所以想通过数值解检验一下解析解的正确性
发表于 2011-1-8 20:24 | 显示全部楼层
[t,y]=ode45(@bcs,[0:1:1],y0);不知道什么意思,步长定值啊?
这个tspan应该改为[0:1:10]吗,
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-11 08:18 , Processed in 0.080003 second(s), 22 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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