运动微分方程中的积分项系数（和t有关）应该如何处理
运动微分方程如图所示其中力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有关，这样的方程变化系数需要通过什么样的处理来求解方程！ 因为积分里有变量t，用quadl求解，太复杂的式子就不行了。看你的K系数并不复杂，试试先手工积分出来，然后再代入方程中。 谢谢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=;
period=2*pi/omega;
=ode45(@diff_equ,,y0)
这部分也是无法得出解的。
所以，对于手工积分系数K虽然能够进行，但是不明白怎样同function以及ode45结合起来。
我这两天看了一些帖子，上面说变系数的微分方程似乎用ode45无法求解。还有的想法是采用内嵌函数，我还是有一些糊涂，主要想知道能够用ode45求解；如果可行的话，应该怎样操作。
希望messenger和大家都能帮忙看一下，毕竟方程推导出来无法求解是一件相当郁闷的事情！ 可以利用嵌套函数共享参数，结构如下，楼主先按我给你的框架自己尝试下，有什么问题再发上来讨论。
function chunshui2003
%这里把Rg,L1，Lambda等一系列已知参数赋好值
function uu=diff_equ(t,u)
uu=zeros(4,1);
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);
%K12,K_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=;
period=2*pi/omega;
=ode45(@diff_equ,,y0);
end
吴老师，你好，根据你给的提示，建立了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=;
period=2*pi/13.1;
=ode45(@diff_equ,,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
=ode45(@diff_equ,,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)中，同样会出现类似的提示。
请你看一下！ 另外，根据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，只是起止时刻为一个步长间隔，将每一次计算后的结果作为下一个步长的初值，并保存每一次计算的最终结果，如此循环直到最后一个时刻完结为止。
按照这种方法微分方程可以计算，请问老师你觉得这个可行吗？ 您那样构造匿名函数不可以的，按下面代码就可以了，注意看下和你原来代码的区别，不会报错，但是有warning，(K12=quadgk(D12,0,2*pi);和后面的几句)估计和你的参数有关，我没有仔细看，你设置断点看看。看看方程的参数范围对求出的积分影响。(上传代码的时候代码放在【code】 【/code】之间(上述【】为[])，这样可读性好的多，还有等号两边预留空格增加程序可读性)
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 = @(alpha) Fsm.*cos(omega_f.*t-p.*alpha); %积分大括号内第一项
C = @(alpha) Fjm.*cos(omega_f.*t-p.*alpha+theta+phi+pi/2); %积分大括号内第二项
C11 = @(alpha) 1+cos(2.*alpha);
C12 = @(alpha) sin(2.*alpha);
C22 = @(alpha) 1-cos(2.*alpha);
C_1 = @(alpha) cos(alpha);
C_2 = @(alpha) sin(alpha); %不同电磁刚度项的积分内容
%
D11 = @(alpha) coeff_1.*C11(alpha).*(B(alpha)+C(alpha)).^2;
D12 = @(alpha) coeff_1.*C12(alpha).*(B(alpha)+C(alpha)).^2;
D22 = @(alpha) coeff_1.*C22(alpha).*(B(alpha)+C(alpha)).^2;
D_1 = @(alpha) coeff_2.*C_1(alpha).*(B(alpha)+C(alpha)).^2;
D_2 = @(alpha) coeff_2.*C_2(alpha).*(B(alpha)+C(alpha)).^2; %电磁刚度积分表达式
%
K11 = quadgk(D11,0,2*pi);
K12 = quadgk(D12,0,2*pi);
K21 = quadgk(D12,0,2*pi);
K22 = quadgk(D22,0,2*pi);
K_1 = quadgk(D_1,0,2*pi);
K_2 = quadgk(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 = ;
period = 2*pi/13.1;
= ode15s(@diff_equ,,y0);
end
[ 本帖最后由 rocwoods 于 2010-7-17 21:55 编辑 ] 吴老师，谢谢你的建议，关于“=”和代码的问题都是我以前没有注意到的，以后会留心的。
程序我运行了，虽然出现warning提示，但可以求解。不过一个问题是，在运行“chunshui2003”这个文件后，计算结果却没有显示。如果是常微分方程组，那么求解命令 = ode15s(@diff_equ,,y0); 一般是在主窗口中运行，但现在是在函数中，结果无法显示。
因为内嵌了子函数“ function uu = diff_equ(t,u) ”，并且主函数没有声明变量，所以对于输出结果这一部分有一些糊涂。
请问老师，“chunshui2003”这个函数的输出结果时间 t 和序列 u 怎样才能够在workspace中出现。 老师，按照你的建议在代码前后分别加了 和，但是显示的结果不是想象中的那样，请老师指正一下。 在chunshui2003前面加个 = 就行了，如下:function = chunshui2003按这样格式：你的代码[反斜杠code]
出现warning说明计算的结果可能不可信，你还是先设置断点自行检查参数的问题吧 恩，谢谢，吴老师，麻烦你了！
页:
[1]