高手救急@求解变参数常微分方程

ljunw003

问题如下，在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

ljunw003

这种方程用ODE45能解吗，如果不能的话又用什么方法呢？
期待高手指教

ChaChing

建議樓主看下本版規則!

ljunw003

楼上说的在理，正在学习论坛版规，虚心向各位请教:handshake

sogooda

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

复制代码

ljunw003

谢谢四楼sogooda朋友的热心帮助。正在学习你的方法，我课题里面出现的变参数有十多个，有问题再向你请教，谢谢你，你让我看到了解决这个问题的希望。

ljunw003

刚才我计算了一下这个方程的解析解，发现跟四楼的方法得到的数值解相差很大，计算数值解时我把xdot222放到m文件里面，
输入的命令是[T,Y] = ode45(@xdot222,[0:10],[1 2])
两个计算结果附在下面了

自变量	参数			解析解		数值解
t	a	b	c	x	y	x	y
0	1	2	1	1	2	0	0
1	1.1	6	4	13.9359807	27.221748	10	10
2	1.2	3	2	36.08111733	50.46852889	40	70
3	1.3	5	2	235.5657149	339.7934075	140	230
4	1.4	3	1	847.7687091	598.2633636	560	660
5	1.5	7	3	10240.90702	20411.81403	2390	3460
6	1.6	2	1	33219.50852	20756.06783	11290	13440
7	1.7	3	2	407147.1404	478972.7534	58810	60110
8	1.8	5	4	6777613.03	15061273.18	338500	578300
9	1.9	4	2	82896088.81	87259003.9	2152660	3211940
10	2	2	3	970330389.8	1455495555	15130610	20734240

期待高手解答

sogooda

先看看趋势对不对，如果只是精度不好，可以通过odeset来改善。
doc odeset

ljunw003

趋势是对的，跟解析解有相似的趋势
应该是自己水平太菜的缘故，按照sogooda朋友提供了思路我修改了半天精度还是没有大的改善，期待高手指点一下options怎样设置才使数值解跟解析解最接近。

ljunw003

感谢大家对我的提问的关注，特别要感谢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，遗憾

weiniuzhu

高手，我也遇到这样一个问题，是关于曲轴连杆动力学方程的，因为压力p也是个变参数，所以想通过数值解检验一下解析解的正确性

weiniuzhu

[t,y]=ode45(@bcs,[0:1:1],y0);不知道什么意思，步长定值啊？
这个tspan应该改为[0:1:10]吗，