运动微分方程中的积分项系数（和t有关）应该如何处理

chunshui2003

运动微分方程如图所示




其中力Fx,Fy,
电磁刚度系数K11,K12,K21,K22,K_1,K_2的表达式如图所示。
其余的参数均为已知量。
如果电磁刚度K系列为常量，那么这是一个简单的运动微分方程，用Rugga—Kutta法可以很容易求解。但是方程中的系数出现了积分项，电磁刚度是随着时间而变化的。加入转频omega=10Hz，则运动周期为1/10s，在一个周期内分成100点，一共求解100个周期，那么电磁刚度K在每一个分解点中都要求解。


我的问题是：应该采用什么样的办法计算电磁刚度K，从而使得运动微分方程能够通过Ode45来求解。
我个人的想法是在设置运动微分方程函数function 时，在里面嵌入quadl命令，求解每一个分解点的电磁刚度，不知道是否可行。请大家帮忙看一下。

也许我说的有一些繁琐，简单来说就是运动微分方程中有的系数不是常量，需要通过积分求解，并且和时间t有关，这样的方程变化系数需要通过什么样的处理来求解方程！

messenger

因为积分里有变量t，用quadl求解，太复杂的式子就不行了。看你的K系数并不复杂，试试先手工积分出来，然后再代入方程中。

chunshui2003

谢谢messenger的意见。求解K系数的积分是没有问题的，只不过t是变化的，采用ode45计算的话，其运动方程中的系数都要为已知。如果在function函数中将运动微分方程写入：
function uu=diff_equ(t,u)

uu=zeros(4,1);
uu(1)=u(2);
uu(2)=-De/(mr*omega)*u(2)-(K11+Ke)/(mr*omega^2)*u(1)-K12/(mr*omega^2)*u(3)+Fx/(mr*delta_0*omega^2)-K_1/(mr*delta_0*omega^2);
uu(3)=u(4);
uu(4)=-De/(mr*omega)*u(4)-K21/(mr*omega^2)*u(1)-(K22+Ke)/(mr*omega^2)*u(3)+Fy/(mr*delta_0*omega^2)-K_2/(mr*delta_0*omega^2);

%K11=@(alpha)quadl(Rg*L1*Lambda_0/(2*sigma^2).*(1+cos(2*alpha)).*(Fsm*cos(omega_f*t-p*alpha)+Fjm*cos(omega_f*t-p*alpha+theta+phi+pi/2).^2)

在这里面我不清楚该如何写入系数K的积分表达式，因为在function里t是不能赋值的，那么系数K就没办法得到具体的数值

在主窗口中：
y0=[0.01 0.01 0.01 0.01];
period=2*pi/omega;
[t,u]=ode45(@diff_equ,[0:period/100:100*period],y0)

这部分也是无法得出解的。

所以，对于手工积分系数K虽然能够进行，但是不明白怎样同function以及ode45结合起来。
我这两天看了一些帖子，上面说变系数的微分方程似乎用ode45无法求解。还有的想法是采用内嵌函数，我还是有一些糊涂，主要想知道能够用ode45求解；如果可行的话，应该怎样操作。

希望messenger和大家都能帮忙看一下，毕竟方程推导出来无法求解是一件相当郁闷的事情！

rocwoods

可以利用嵌套函数共享参数，结构如下，楼主先按我给你的框架自己尝试下，有什么问题再发上来讨论。

1. 
   
2. function chunshui2003
   
3. %这里把Rg,L1，Lambda等一系列已知参数赋好值
   
4. function uu=diff_equ(t,u)
   
5. uu=zeros(4,1);
   
6. K11=quadl(@(alpha) Rg*L1*Lambda_0/(2*sigma^2).*(1+cos(2*alpha)).*(Fsm*cos(omega_f*t-p*alpha)+Fjm*cos(omega_f*t-p*alpha+theta+phi+pi/2).^2),0,2*pi);
   
7. %K12,K_1等等类似写出来
   
8. uu(1)=u(2);
   
9. uu(2)=-De/(mr*omega)*u(2)-(K11+Ke)/(mr*omega^2)*u(1)-K12/(mr*omega^2)*u(3)+Fx/(mr*delta_0*omega^2)-K_1/(mr*delta_0*omega^2);
   
10. uu(3)=u(4);
    
11. uu(4)=-De/(mr*omega)*u(4)-K21/(mr*omega^2)*u(1)-(K22+Ke)/(mr*omega^2)*u(3)+Fy/(mr*delta_0*omega^2)-K_2/(mr*delta_0*omega^2);
    
12. end
    
13. y0=[0.01 0.01 0.01 0.01];
    
14. period=2*pi/omega;
    
15. [t,u]=ode45(@diff_equ,[0:period/100:100*period],y0);
    
16. 
    
17. end
    
18. 
    

复制代码

chunshui2003

吴老师，你好，根据你给的提示，建立了chunshui2003.m文件如下：


function chunshui2003

%%%% 基本参数数值 %%%%%
mr=600*10^3;            
De=0.5*10^6;            
Ke=0.5*10^10;           
delta_0=18*10^-3;      
e_0=2*10^-3;            
omega=13.1;              %水轮发电机额定转速，单位rad/s


%%%% 电磁刚度参数数值 %%%%%
Rg=6.24;      
L1=2.1;      
k_u=1.102;   
mu_0=4*pi*10^(-7);   
f=50;        
theta=30.64/180*pi;
phi=acos(0.875);   
p=40;               
Fsm=19210;         
Fjm=24214;         
Lambda_0=mu_0./(k_u*delta_0);      
sigma=k_u*delta_0;      
omega_f=2*pi*f./p;      


function uu=diff_equ(t,u)

%


%%%%% Fx,Fy 表达式 %%%%%%

Fx=mr.*omega.^2.*e_0.*cos(omega.*t);
Fy=mr.*omega.^2.*e_0.*sin(omega.*t);


%%%% 电磁刚度表达式 %%%%%

coeff_1=Rg.*L1.*Lambda_0./(2*sigma.^2);
coeff_2=coeff_1.*sigma;

B=Fsm.*cos(omega_f.*t-p.*alpha);                     %积分大括号内第一项
C=Fjm.*cos(omega_f.*t-p.*alpha+theta+phi+pi/2);      %积分大括号内第二项

C11=1+cos(2.*alpha);
C12=sin(2.*alpha);
C22=1-cos(2.*alpha);
C_1=cos(alpha);
C_2=sin(alpha);                  %不同电磁刚度项的积分内容
%
D11=coeff_1.*C11.*(B+C).^2;
D12=coeff_1.*C12.*(B+C).^2;
D22=coeff_1.*C22.*(B+C).^2;
D_1=coeff_2.*C_1.*(B+C).^2;
D_2=coeff_2.*C_2.*(B+C).^2;                   %电磁刚度积分表达式
%
K11=quadl(@(alpha)D11,0,2*pi);
K12=quadl(@(alpha)D12,0,2*pi);
K21=quadl(@(alpha)D12,0,2*pi);
K22=quadl(@(alpha)D22,0,2*pi);
K_1=quadl(@(alpha)D_1,0,2*pi);
K_2=quadl(@(alpha)D_2,0,2*pi);              %电磁刚度求解



uu=zeros(4,1);
uu(1)=u(2);
uu(2)=-De/(mr*omega)*u(2)-(K11+Ke)/(mr*omega^2)*u(1)-K12/(mr*omega^2)*u(3)+Fx/(mr*delta_0*omega^2)-K_1/(mr*delta_0*omega^2);
uu(3)=u(4);
uu(4)=-De/(mr*omega)*u(4)-K21/(mr*omega^2)*u(1)-(K22+Ke)/(mr*omega^2)*u(3)+Fy/(mr*delta_0*omega^2)-K_2/(mr*delta_0*omega^2);
end


y0=[0.001 0.001 0.001 0.001];
period=2*pi/13.1;
[t,u]=ode45(@diff_equ,[0:period/100:100*period],y0);
end

运行后在主窗口出现错误提示：
??? Error using ==> alpha
Too many output arguments.

Error in ==> chunshui2003>diff_equ at 45
B=Fsm.*cos(omega_f.*t-p.*alpha);                     %积分大括号内第一项

Error in ==> odearguments at 110
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in ==> chunshui2003 at 77
[t,u]=ode45(@diff_equ,[0:period/100:100*period],y0);


如果将参数部分
%%%% 基本参数数值 %%%%%
mr=600*10^3;            
De=0.5*10^6;            
Ke=0.5*10^10;           
delta_0=18*10^-3;      
e_0=2*10^-3;            
omega=13.1;              %水轮发电机额定转速，单位rad/s


%%%% 电磁刚度参数数值 %%%%%
Rg=6.24;      
L1=2.1;      
k_u=1.102;   
mu_0=4*pi*10^(-7);   
f=50;        
theta=30.64/180*pi;
phi=acos(0.875);   
p=40;               
Fsm=19210;         
Fjm=24214;         
Lambda_0=mu_0./(k_u*delta_0);      
sigma=k_u*delta_0;      
omega_f=2*pi*f./p;      

写入到子函数function uu=diff_equ(t,u)中，同样会出现类似的提示。
请你看一下！

chunshui2003

另外，根据messenger的提示，我昨天想到的一个方法是，对电磁刚度系数K系列先采用积分，比如K11：
%%%%%不同时刻t下电磁刚度数值K11的计算  %%%%%

function y=jifen_k11(t)              
function k11=k11(alpha)              
Rg=6.24;      
L_1=2.1;      
k_u=1.102;   
mu_0=4*pi*10^(-7);     
delta_0=18e-003;      
f=50;         
theta=30.64/180*pi;
phi=acos(0.875);   
p=40;               
Fsm=19210;         
Fjm=24214;         
Lambda_0=mu_0./(k_u*delta_0);      
sigma=k_u*delta_0;         
omega_f=2*pi*f./p;        %发电机同步转速

Coeff_1=Rg.*L_1.*Lambda_0./(2.*sigma.^2);    %%%%第一类电磁刚度前面的常系数

coeff_k11=1+cos(2.*alpha);             %不同的内在表达式
A_1=Fsm*cos(omega_f.*t-p.*alpha);
A_2=Fjm*cos(omega_f.*t-p.*alpha+theta+phi+pi./2);
k11=Coeff_1.*coeff_k11.*(A_1+A_2).^2;            %K11的积分表达式
end
y=quadl(@k11,0,2*pi);
end

主窗口输入：
omega=13.1;               
period=2*pi/omega;        
tt=(0:period/100:100*period);            
y=zeros(size(tt));                       
for ii=1:length(tt)
    tt(ii)
    y(ii)=jifen_k11(tt(ii));
end

其余的刚度系数类似，然后将它们赋给总体刚度系数K，形成一个矩阵以备调用。
在进行微分方程计算时，将每个电磁刚度系数都作为全局变量，在不同的t时刻由总体刚度矩阵K赋给不同的值，仍然采用ode45，只是起止时刻为一个步长间隔，将每一次计算后的结果作为下一个步长的初值，并保存每一次计算的最终结果，如此循环直到最后一个时刻完结为止。

按照这种方法微分方程可以计算，请问老师你觉得这个可行吗？

rocwoods

您那样构造匿名函数不可以的，按下面代码就可以了，注意看下和你原来代码的区别，不会报错，但是有warning，(K12=quadgk(D12,0,2*pi);和后面的几句)估计和你的参数有关，我没有仔细看，你设置断点看看。看看方程的参数范围对求出的积分影响。(上传代码的时候代码放在【code】 【/code】之间(上述【】为[])，这样可读性好的多，还有等号两边预留空格增加程序可读性)

1. 
   
2. function chunshui2003
   
3. 
   
4. %%%% 基本参数数值 %%%%%
   
5. mr = 600*10^3;            
   
6. De = 0.5*10^6;            
   
7. Ke = 0.5*10^10;           
   
8. delta_0 = 18*10^-3;      
   
9. e_0 = 2*10^-3;            
   
10. omega = 13.1;              %水轮发电机额定转速，单位rad/s
    
11. 
    
12. 
    
13. %%%% 电磁刚度参数数值 %%%%%
    
14. Rg = 6.24;      
    
15. L1 = 2.1;      
    
16. k_u = 1.102;   
    
17. mu_0 = 4*pi*10^(-7);   
    
18. f = 50;        
    
19. theta = 30.64/180*pi;
    
20. phi = acos(0.875);   
    
21. p = 40;               
    
22. Fsm = 19210;         
    
23. Fjm = 24214;         
    
24. Lambda_0 = mu_0./(k_u*delta_0);      
    
25. sigma = k_u*delta_0;      
    
26. omega_f = 2*pi*f./p;      
    
27. 
    
28. 
    
29. function uu = diff_equ(t,u)
    
30. %%%%% Fx,Fy 表达式 %%%%%%
    
31. 
    
32. Fx = mr.*omega.^2.*e_0.*cos(omega.*t);
    
33. Fy = mr.*omega.^2.*e_0.*sin(omega.*t);
    
34. 
    
35. 
    
36. %%%% 电磁刚度表达式 %%%%%
    
37. 
    
38. coeff_1 = Rg.*L1.*Lambda_0./(2*sigma.^2);
    
39. coeff_2 = coeff_1.*sigma;
    
40. 
    
41. B = @(alpha) Fsm.*cos(omega_f.*t-p.*alpha);                     %积分大括号内第一项
    
42. C = @(alpha) Fjm.*cos(omega_f.*t-p.*alpha+theta+phi+pi/2);      %积分大括号内第二项
    
43. 
    
44. C11 = @(alpha) 1+cos(2.*alpha);
    
45. C12 = @(alpha) sin(2.*alpha);
    
46. C22 = @(alpha) 1-cos(2.*alpha);
    
47. C_1 = @(alpha) cos(alpha);
    
48. C_2 = @(alpha) sin(alpha);                  %不同电磁刚度项的积分内容
    
49. %
    
50. D11 = @(alpha) coeff_1.*C11(alpha).*(B(alpha)+C(alpha)).^2;
    
51. D12 = @(alpha) coeff_1.*C12(alpha).*(B(alpha)+C(alpha)).^2;
    
52. D22 = @(alpha) coeff_1.*C22(alpha).*(B(alpha)+C(alpha)).^2;
    
53. D_1 = @(alpha) coeff_2.*C_1(alpha).*(B(alpha)+C(alpha)).^2;
    
54. D_2 = @(alpha) coeff_2.*C_2(alpha).*(B(alpha)+C(alpha)).^2;                   %电磁刚度积分表达式
    
55. %
    
56. K11 = quadgk(D11,0,2*pi);
    
57. K12 = quadgk(D12,0,2*pi);
    
58. K21 = quadgk(D12,0,2*pi);
    
59. K22 = quadgk(D22,0,2*pi);
    
60. K_1 = quadgk(D_1,0,2*pi);
    
61. K_2 = quadgk(D_2,0,2*pi);              %电磁刚度求解
    
62. 
    
63. uu = zeros(4,1);
    
64. uu(1) = u(2);
    
65. uu(2) = -De/(mr*omega)*u(2)-(K11+Ke)/(mr*omega^2)*u(1)-K12/(mr*omega^2)*u(3)+Fx/(mr*delta_0*omega^2)-K_1/(mr*delta_0*omega^2);
    
66. uu(3) = u(4);
    
67. uu(4) = -De/(mr*omega)*u(4)-K21/(mr*omega^2)*u(1)-(K22+Ke)/(mr*omega^2)*u(3)+Fy/(mr*delta_0*omega^2)-K_2/(mr*delta_0*omega^2);
    
68. end
    
69. 
    
70. 
    
71. y0 = [0.001 0.001 0.001 0.001];
    
72. period = 2*pi/13.1;
    
73. [t,u] = ode15s(@diff_equ,[0:period/100:100*period],y0);
    
74. end
    
75. 
    
76. 
    

复制代码

[ 本帖最后由 rocwoods 于 2010-7-17 21:55 编辑 ]

chunshui2003

吴老师，谢谢你的建议，关于“=”和代码的问题都是我以前没有注意到的，以后会留心的。
程序我运行了，虽然出现warning提示，但可以求解。不过一个问题是，在运行“chunshui2003”这个文件后，计算结果却没有显示。如果是常微分方程组，那么求解命令[code ] [t,u] = ode15s(@diff_equ,[0:period/100:100*period],y0); [/code]一般是在主窗口中运行，但现在是在函数中，结果无法显示。
因为内嵌了子函数“ function uu = diff_equ(t,u) ”，并且主函数没有声明变量，所以对于输出结果这一部分有一些糊涂。
请问老师，“chunshui2003”这个函数的输出结果时间 t 和序列 u 怎样才能够在workspace中出现。

chunshui2003

老师，按照你的建议在代码前后分别加了

1. 和

复制代码

，但是显示的结果不是想象中的那样，请老师指正一下。

rocwoods

在chunshui2003前面加个[t,u] = 就行了，如下:

1. function [t,u] = chunshui2003

复制代码

按这样格式：[code]你的代码[反斜杠code]
出现warning说明计算的结果可能不可信，你还是先设置断点自行检查参数的问题吧

chunshui2003

恩，谢谢，吴老师，麻烦你了！