matlab中参数的敏感性分析

chunshui2003

function uu=yundongguiji(t,u)

m1=1.5*10^4; m2=1.1*10^4; c1=6.0*10^4; c2=7.0*10^4;            
k1=8.5*10^7; k2=6.5*10^7; k3=3.5*10^7;         
delta0=8*10^-3; e2=0.3*10^-3; mu0=4*pi*10^-7;            
kj=5.2; R=1.2; L=0.5; omega=10;               

K1=25*k1+9*k2+k3; K2=-5*k1+3*k2+3*k3; K3=k1+k2+9*k3;

e1=0.5*10-3;           %发电机转子偏心距离

B1=e1*omega^2*cos(omega*t)/delta0;
B2=e1*omega^2*sin(omega*t)/delta0;
B3=e2*omega^2*cos(omega*t)/delta0;
B4=e2*omega^2*sin(omega*t)/delta0;

Ij=1065;               %发电机转子励磁电流

F0=R*L*pi*kj^2*Ij^2*mu0/(m1*delta0^3);

Z1=(1-sqrt(1-u(1).^2-u(3).^2))/sqrt(u(1).^2+u(3).^2);
Z2=sqrt(u(5).^2+u(7).^2)/sqrt(u(1).^2+u(3).^2);

uu=zeros(8,1);
uu(1)=u(2);
uu(2)=B1+F0./(1-u(1).^2-u(3).^2)*(Z1+Z1.^3+Z1.^5).*u(1)./sqrt(u(1).^2+u(3).^2)-c1./m1.*u(2)-u(1).*(K1+K2.*Z2)./(16.*m1);
uu(3)=u(4);
uu(4)=B2+F0./(1-u(1).^2-u(3).^2)*(Z1+Z1.^3+Z1.^5).*u(3)./sqrt(u(1).^2+u(3).^2)-c1./m1.*u(4)-u(3).*(K1+K2.*Z2)./(16.*m1);
uu(5)=u(6);
uu(6)=B3-c2./m2.*u(6)-u(5).*(K3+K2.*Z2)./(16.*m2);
uu(7)=u(8);
uu(8)=B4-c2./m2.*u(8)-u(7)*(K3+K2.*Z2)./(16.*m2);

程序可能有点长，但关键的就是最后面的几行，剩下的基本都是参数。
用ode45命令可以求出相应的u(1)-u(8)的数值，是关于时间t的。
但是我希望得到的是u(1)-u(8)关于Ij的敏感性关系，也就是Ij从1050变化到1120的过程中，u对应于Ij的曲线。
有一种笨方法是每次给出一个Ij的值，然后求解u值，但这样太麻烦了。
请问，如何可以在程序中实现随着Ij的变化，得到相应的u值。
希望能够赋给程序，因为我现在确实无法想出来，谢谢。

[ 本帖最后由 ChaChing 于 2009-12-1 19:57 编辑 ]

ChaChing

LZ解决否? 看看这个是否有用
关于求解变参数微分方程
http://forum.vibunion.com/forum/ ... id=44972&page=1

chunshui2003

感谢ChaChing，原来我的思路有一些问题，现在基本解决了。

ChaChing

建议与大家说明下, 分享你的成果, 做个结束!

chunshui2003

之前思路上存在一些问题。微分方程组某一个参数的变化有可能导致系统的特性改变，由于采用Rugga-Kutta法，位移、速度及加速度等均是和时间有关的函数，不同的时刻他们的数值是发生变化的。如果希望了解方程组中某参数的变化对整体系统特性的影响，那么只能知道某一个具体时刻的情况，如t=0.001s或者t=0.5s等等。所以当我看到某些文献上做参数敏感性分析得到是由单一点组成的直线或曲线时也想尝试一下，但越发觉得不对，应该说是没有考虑到具体时间这一因素。而采用分岔图来说明就比较合理了，毕竟是周期性取点。
所以最终我选择了参数在一定范围内以某间隔变化，观察系统的轨迹、时频域曲线等，而非分岔图和poincare映射图。在此我要说明一下，论坛里有好多朋友是搞汽轮机和航天发动机转子方面的，这类转子的转速较高，容易发生分岔现象，而我是搞水轮机的，转速较低，上述现象较难发生，但也不是一定不会发生，所以侧重于其他方面更好一些。
说了很多，关键就是想求出系统随参数变化的由单一点组成的线性或类线性曲线是不可能的，因为时间是变化的。
也许我的理解存在一些问题，欢迎大家多提意见。