您所在位置：主题：[讨论]Fortran转换matlab计算微分方程组，提示：可能为刚性问题
chunshui2003
[专家分：0] 发布于 2010-10-31 11:01:00
matlab的代码如下：
[code]
function uu=differential_equation(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; %机组旋转频率 单位Hz
K1=25*k1+9*k2+k3;
K2=-5*k1+3*k2+3*k3;
K3=k1+k2+9*k3; %% 自定义的系数，方便使用
e1=0.2*10^-3; %发电机转子偏心距离
B1=e1*omega^2*cos(2*pi*omega*t)/delta0;
B2=e1*omega^2*sin(2*pi*omega*t)/delta0;
Ij=500; %% 励磁电流
%%%%%%%%%%%%%% 线性密封力参数 %%%%%%%%%%%%%%%%%%
alpha0=1.0; % 动量修正系数
rho=1.0; % 液体密度，单位kg/m^3
lambda=1.0; % 沿程损失系数
xi=1.5; % 局部损失系数
Rm=2.925; % 密封半径 单位m
lm=0.43; % 密封长度 单位m
v=3.537; % 液体轴向流速 单位m/s
cm=2.5*10^-3; % 密封间隙 单位m
Kxx=lambda*lm.^2*rho*pi*v.^2*Rm*(1+xi)/(8*cm.^2*(1+xi+lambda*lm/2/cm));
Kyy=Kxx;
Kyx=-(Rm*omega/2)*(pi*rho*v*lm.^2/2/cm)*(1+xi+lambda*lm/6/cm);
Kxy=-Kyx; %%% 刚度系数
xishu=(2*(1+xi)+lambda*lm/4/cm)/(3*(1+xi)+lambda*lm/2/cm);
Dxx=(rho*v*lm.^2/2/cm)*pi*Rm*(1+xi+lambda*lm/6/cm);
Dyy=Dxx;
Dyx=-(Rm*omega/2)*(pi*rho*alpha0*lm.^3/cm)*xishu;
Dxy=-Dyx; %%%%% 阻尼系数
mxx=pi*Rm*rho*alpha0*lm.^3/cm*xishu;
myy=mxx; %%%% 惯性系数
%%%%%%%%%% 线性密封力合成参数 %%%%%%%%%%%%
M2x=m2+mxx;
M2y=m2+myy;
D1=(c2+Dxx)/M2x;
D2=Dxy/M2x;
D3=(c2+Dyy)/M2y;
D4=Dyx/M2y;
Kg1=Kxy/M2x;
Kg2=Kyx/M2y;
B3=m2*e2*omega.^2*cos(2*pi*omega*t)/(M2x*delta0);
B4=m2*e2*omega.^2*sin(2*pi*omega*t)/(M2y*delta0);
%%%%%%%%%%% 涉及到变量的一些参数
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);
Z3=1/Z2;
F0=R*L*pi*kj.^2*Ij.^2*mu0/(m1*delta0.^3);
F1=(Z1+Z1^3+Z1^5)/(1-u(1)^2-u(3)^2);
F=F0*F1;
%%%%%%%%% 转子-转轮系统运动微分方程
uu=zeros(8,1);
uu(1)=u(2);
uu(2)=B1+F*u(1)/sqrt(u(1)^2+u(3)^2)-c1/m1*u(2)-1/(16*m1)*(K1+K2*Z2)*u(1);
uu(3)=u(4);
uu(4)=B2+F*u(3)/sqrt(u(1)^2+u(3)^2)-c1/m1*u(4)-1/(16*m1)*(K1+K2*Z2)*u(3);
uu(5)=u(6);
uu(6)=B3-D1*u(6)-D2*u(8)-1/(16*M2x)*(K3+K2*Z3)*u(5)-Kg1*u(7);
uu(7)=u(8);
uu(8)=B4-D3*u(8)-D4*u(6)-1/(16*M2y)*(K3+K2*Z3)*u(7)-Kg2*u(5);
[/code]
求解：
[code]
clc
clear
tic
y0=[0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 ];
period=1/10; %采样的周期，转速的单位是10Hz,所以周期为“1/10”；如果转速的单位为“rad/s”，那么周期为“2*pi/omega”
[t,u]=ode45(@differential_equation,[0:period/100:100*period],y0);
toc
[/code]
大概计算时间为2.5秒，没有问题，非刚性。
改为Fortran程序如下：
[code]
program main
use IMSL
implicit none
integer MXPARM , N
parameter (MXPARM = 50, N = 8)
integer, parameter :: MABSE = 1, MBDF = 2, MSLOVE = 2 ! 选择绝对误差，采用BDF算法，Jacobian矩阵由机器提供
integer IDO, ISTEP, NOUT
real A(1,1), PARAM(MXPARM), T, TEND, TOL, Y(N)
external FCN, FCNJ
T = 0.0
Y(1) = 1.0E-4
Y(2) = 1.0E-4
Y(3) = 1.0E-4
Y(4) = 1.0E-4
Y(5) = 1.0E-4
Y(6) = 1.0E-4
Y(7) = 1.0E-4
Y(8) = 1.0E-4
TOL = 0.001 ! 绝对误差
PARAM = 0 ! 采用默认值
PARAM(4) = 1000000 ! 允许的最大步数
PARAM(10) = MABSE
PARAM(12) = MBDF
PARAM(13) = MSLOVE ! 1是用户提供Jacobian，2是差分Jacobian
write(NOUT,99998)
IDO = 1
do ISTEP = 1,100
TEND = ISTEP
call IVPAG (IDO, N, FCN, FCNJ, A, T, TEND, TOL, PARAM, Y)
write(NOUT, '(I6, 9F12.3)') ISTEP, T, Y
end do
call IVPAG (3, N, FCN, FCNJ, A, T, TEND, TOL, PARAM, Y) ! 释放内存
99998 FORMAT (11X, 'T', 14X, 'Y(1)', 11X, 'Y(2)', 11X, 'Y(3)', 11X, 'Y(4)', 11X, 'Y(5)', 11X, 'Y(6)', 11X, 'Y(7)', 11X, 'Y(8)')
stop
end program
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine FCN (N, T, Y, YPRIME)
integer N
real T, Y(N), YPRIME(N)
REAL, PARAMETER :: PI = 3.14159
real F0, F1, F, e
save F0
real, parameter :: Ij = 500.0
real K_1, K_2, K_3
save K_1, K_2, K_3
real B1, B2, B3, B4
real Z1, Z2, Z3
real Kxx, Kyy, Kxy, Kyx, xishu, Dxx, Dyy, Dxy, Dyx, mxx, myy
save Kxx, Kyy, Kxy, Kyx, xishu, Dxx, Dyy, Dxy, Dyx, mxx, myy
real M2x, M2y, D1, D2, D3, D4, Kg1, Kg2
save M2x, M2y, D1, D2, D3, D4, Kg1, Kg2
real m1, m2, c1, c2, k1, k2, k3, delta0, e1, e2, mu0, kj, R, L, omega !运动微分方程参数
save m1, m2, c1, c2, k1, k2, k3, delta0, e1, e2, mu0, kj, R, L, omega !使其成为全局常数，保持不变
real alpha0, rho, lambda, xi, Rm, lm, v, cm !密封参数
save alpha0, rho, lambda, xi, Rm, lm, v, cm
!!!!! 运动微分方程参数 !!!!!!
m1 = 1.5E4 !发电机转子质量
m2 = 1.1E4 !水轮机转轮质量
c1 = 6.0E4; !发电机转子阻尼系数
c2 = 7.0E4; !水轮机转轮阻尼系数
k1 = 8.5E7; !上导轴承刚度系数
k2 = 6.5E7; !下导轴承刚度系数
k3 = 3.5E7; !水导轴承刚度系数
delta0 = 8.0E-3; !发电机转子不偏心时平均气隙长度
e1 = 0.2E-3 !发电机转子偏心距离
e2 = 0.3E-3; !水轮机转轮偏心距离
mu0 = 4.0*PIE-7; !空气磁导系数
kj = 5.2E0; !气隙基波磁动势系数
R = 1.2E0; !发电机转子半径
L = 0.5E0; !发电机转子长度
omega = 10.0E0; !机组旋转频率 单位是Hz
!!!!! 转轮密封力参数 !!!!!!
alpha0 = 1.0E0; !动量修正系数
rho = 1.0E0; !液体密度，单位kg/m^3
lambda = 1.0E0; ! 沿程损失系数
xi = 1.5E0; ! 局部损失系数
Rm = 2.925E0; ! 密封半径 单位m
lm = 0.43E0; ! 密封长度 单位m
v = 3.537E0; ! 液体轴向流速 单位m/s
cm = 2.5E-3; ! 密封间隙 单位m
K_1 = 25*k1 + 9*k2 + k3
K_2 = -5*k1 + 3*k2 + 3*k3
K_3 = k1 + k2 + 9*k3
B1 = e1*omega**2 * cos(2*PI*omega*T)/delta0
B2 = e1*omega**2 * sin(2*PI*omega*T)/delta0
Kxx = lambda*lm**2 *rho*PI*v**2 *Rm*(1.0+xi)/(8.0*cm**2 *(1.0+xi+lambda*lm/2.0/cm))
Kyy = Kxx
Kyx = -(Rm*omega/2.0)*(PI*rho*v*lm**2 /2.0/cm)*(1.0+xi+lambda*lm/6.0/cm)
Kxy = -Kyx ! 刚度系数
xishu = (2.0*(1+xi)+lambda*lm/4.0/cm)/(3.0*(1.0+xi)+lambda*lm/2.0/cm)
Dxx = (rho*v*lm**2 /2.0/cm)*PI*Rm*(1.0+xi+lambda*lm/6.0/cm)
Dyy = Dxx
Dyx = -(Rm*omega/2.0)*(PI*rho*alpha0*lm**3 /cm)*xishu
Dxy = -Dyx !阻尼系数
mxx = PI*Rm*rho*alpha0*lm**3 /cm*xishu;
myy = mxx ! 惯性系数
M2x = m2 + mxx
M2y = m2 + myy
D1 = (c2 + Dxx)/M2x
D2 = Dxy/M2x
D3 = (c2 + Dyy)/M2y
D4 = Dyx/M2y
Kg1 = Kxy/M2x
Kg2 = Kyx/M2y
B3 = m2*e2*omega**2 * cos(2*pi*omega*T)/(M2x*delta0)
B4 = m2*e2*omega**2 * sin(2*pi*omega*T)/(M2y*delta0)
Z1 = (1.0 - sqrt(1.0 - Y(1)**2 -Y(3)**2 )) / sqrt(Y(1)**2 + Y(3)**2)
Z2 = sqrt(Y(5)**2 + Y(7)**2) / sqrt(Y(1)**2 + Y(3)**2)
Z3 = 1.0 / Z2
!!!!! 不平衡磁拉力 !!!!!!11
!e = sqrt(Y(1)**2 + Y(3)**2)
F0 = R*L*PI*mu0*(Ij*Kj)**2 / (m1*delta0**3)
F1 = (Z1 + Z1**3 + Z1**5) / (1-Y(1)**2 - Y(3)**2)
F = F0 * F1
!!!!!!! 系统运动微分方程 !!!!!!!!!!!!!1
YPRIME(1) = Y(2)
YPRIME(2) = B1 + F*Y(1)/sqrt(Y(1)**2 + Y(3)**2) - c1/m1*Y(2) - 1.0/(16.0*m1) * (K_1+K_2*Z2) * Y(1)
YPRIME(3) = Y(4)
YPRIME(4) = B2 + F*Y(3)/sqrt(Y(1)**2 + Y(3)**2) - c1/m1*Y(4) - 1.0/(16.0*m1) * (K_1+K_2*Z2) * Y(3)
YPRIME(5) = Y(6)
YPRIME(6) = B3 - D1*Y(6) - D2*Y(8) - 1.0/(16.0*M2x) * (K_3+K_2*Z3) * Y(5) - Kg1 * Y(7)
YPRIME(7) = Y(8)
YPRIME(8) = B4 - D3*Y(8) - D2*Y(6) - 1.0/(16.0*M2y) * (K_3+K_2*Z3) * Y(7) - Kg2 * Y(5)
return
end subroutine FCN
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine FCNJ (N, T, Y, DYPDY)
integer N
real T, Y(N), DYPDY(N,*)
return
end subroutine FCNJ ! 实际上并没有用到FCNJ
[/code]
结果提示如图所示：
之前一直用matlab计算微分方程组，所以对fortran不是很了解。但Fortran的计算速度是matlab无法比拟的，所以想修改程序。麻烦大家看一下究竟是什么原因造成这种现象的（方程应该不是刚性的，否则ode45根本不可能在2.5秒之内就计算出来）。
[code]
function uu=differential_equation(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; %机组旋转频率 单位Hz
K1=25*k1+9*k2+k3;
K2=-5*k1+3*k2+3*k3;
K3=k1+k2+9*k3; %% 自定义的系数，方便使用
e1=0.2*10^-3; %发电机转子偏心距离
B1=e1*omega^2*cos(2*pi*omega*t)/delta0;
B2=e1*omega^2*sin(2*pi*omega*t)/delta0;
Ij=500; %% 励磁电流
%%%%%%%%%%%%%% 线性密封力参数 %%%%%%%%%%%%%%%%%%
alpha0=1.0; % 动量修正系数
rho=1.0; % 液体密度，单位kg/m^3
lambda=1.0; % 沿程损失系数
xi=1.5; % 局部损失系数
Rm=2.925; % 密封半径 单位m
lm=0.43; % 密封长度 单位m
v=3.537; % 液体轴向流速 单位m/s
cm=2.5*10^-3; % 密封间隙 单位m
Kxx=lambda*lm.^2*rho*pi*v.^2*Rm*(1+xi)/(8*cm.^2*(1+xi+lambda*lm/2/cm));
Kyy=Kxx;
Kyx=-(Rm*omega/2)*(pi*rho*v*lm.^2/2/cm)*(1+xi+lambda*lm/6/cm);
Kxy=-Kyx; %%% 刚度系数
xishu=(2*(1+xi)+lambda*lm/4/cm)/(3*(1+xi)+lambda*lm/2/cm);
Dxx=(rho*v*lm.^2/2/cm)*pi*Rm*(1+xi+lambda*lm/6/cm);
Dyy=Dxx;
Dyx=-(Rm*omega/2)*(pi*rho*alpha0*lm.^3/cm)*xishu;
Dxy=-Dyx; %%%%% 阻尼系数
mxx=pi*Rm*rho*alpha0*lm.^3/cm*xishu;
myy=mxx; %%%% 惯性系数
%%%%%%%%%% 线性密封力合成参数 %%%%%%%%%%%%
M2x=m2+mxx;
M2y=m2+myy;
D1=(c2+Dxx)/M2x;
D2=Dxy/M2x;
D3=(c2+Dyy)/M2y;
D4=Dyx/M2y;
Kg1=Kxy/M2x;
Kg2=Kyx/M2y;
B3=m2*e2*omega.^2*cos(2*pi*omega*t)/(M2x*delta0);
B4=m2*e2*omega.^2*sin(2*pi*omega*t)/(M2y*delta0);
%%%%%%%%%%% 涉及到变量的一些参数
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);
Z3=1/Z2;
F0=R*L*pi*kj.^2*Ij.^2*mu0/(m1*delta0.^3);
F1=(Z1+Z1^3+Z1^5)/(1-u(1)^2-u(3)^2);
F=F0*F1;
%%%%%%%%% 转子-转轮系统运动微分方程
uu=zeros(8,1);
uu(1)=u(2);
uu(2)=B1+F*u(1)/sqrt(u(1)^2+u(3)^2)-c1/m1*u(2)-1/(16*m1)*(K1+K2*Z2)*u(1);
uu(3)=u(4);
uu(4)=B2+F*u(3)/sqrt(u(1)^2+u(3)^2)-c1/m1*u(4)-1/(16*m1)*(K1+K2*Z2)*u(3);
uu(5)=u(6);
uu(6)=B3-D1*u(6)-D2*u(8)-1/(16*M2x)*(K3+K2*Z3)*u(5)-Kg1*u(7);
uu(7)=u(8);
uu(8)=B4-D3*u(8)-D4*u(6)-1/(16*M2y)*(K3+K2*Z3)*u(7)-Kg2*u(5);
[/code]
求解：
[code]
clc
clear
tic
y0=[0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001 ];
period=1/10; %采样的周期，转速的单位是10Hz,所以周期为“1/10”；如果转速的单位为“rad/s”，那么周期为“2*pi/omega”
[t,u]=ode45(@differential_equation,[0:period/100:100*period],y0);
toc
[/code]
大概计算时间为2.5秒，没有问题，非刚性。
改为Fortran程序如下：
[code]
program main
use IMSL
implicit none
integer MXPARM , N
parameter (MXPARM = 50, N = 8)
integer, parameter :: MABSE = 1, MBDF = 2, MSLOVE = 2 ! 选择绝对误差，采用BDF算法，Jacobian矩阵由机器提供
integer IDO, ISTEP, NOUT
real A(1,1), PARAM(MXPARM), T, TEND, TOL, Y(N)
external FCN, FCNJ
T = 0.0
Y(1) = 1.0E-4
Y(2) = 1.0E-4
Y(3) = 1.0E-4
Y(4) = 1.0E-4
Y(5) = 1.0E-4
Y(6) = 1.0E-4
Y(7) = 1.0E-4
Y(8) = 1.0E-4
TOL = 0.001 ! 绝对误差
PARAM = 0 ! 采用默认值
PARAM(4) = 1000000 ! 允许的最大步数
PARAM(10) = MABSE
PARAM(12) = MBDF
PARAM(13) = MSLOVE ! 1是用户提供Jacobian，2是差分Jacobian
write(NOUT,99998)
IDO = 1
do ISTEP = 1,100
TEND = ISTEP
call IVPAG (IDO, N, FCN, FCNJ, A, T, TEND, TOL, PARAM, Y)
write(NOUT, '(I6, 9F12.3)') ISTEP, T, Y
end do
call IVPAG (3, N, FCN, FCNJ, A, T, TEND, TOL, PARAM, Y) ! 释放内存
99998 FORMAT (11X, 'T', 14X, 'Y(1)', 11X, 'Y(2)', 11X, 'Y(3)', 11X, 'Y(4)', 11X, 'Y(5)', 11X, 'Y(6)', 11X, 'Y(7)', 11X, 'Y(8)')
stop
end program
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine FCN (N, T, Y, YPRIME)
integer N
real T, Y(N), YPRIME(N)
REAL, PARAMETER :: PI = 3.14159
real F0, F1, F, e
save F0
real, parameter :: Ij = 500.0
real K_1, K_2, K_3
save K_1, K_2, K_3
real B1, B2, B3, B4
real Z1, Z2, Z3
real Kxx, Kyy, Kxy, Kyx, xishu, Dxx, Dyy, Dxy, Dyx, mxx, myy
save Kxx, Kyy, Kxy, Kyx, xishu, Dxx, Dyy, Dxy, Dyx, mxx, myy
real M2x, M2y, D1, D2, D3, D4, Kg1, Kg2
save M2x, M2y, D1, D2, D3, D4, Kg1, Kg2
real m1, m2, c1, c2, k1, k2, k3, delta0, e1, e2, mu0, kj, R, L, omega !运动微分方程参数
save m1, m2, c1, c2, k1, k2, k3, delta0, e1, e2, mu0, kj, R, L, omega !使其成为全局常数，保持不变
real alpha0, rho, lambda, xi, Rm, lm, v, cm !密封参数
save alpha0, rho, lambda, xi, Rm, lm, v, cm
!!!!! 运动微分方程参数 !!!!!!
m1 = 1.5E4 !发电机转子质量
m2 = 1.1E4 !水轮机转轮质量
c1 = 6.0E4; !发电机转子阻尼系数
c2 = 7.0E4; !水轮机转轮阻尼系数
k1 = 8.5E7; !上导轴承刚度系数
k2 = 6.5E7; !下导轴承刚度系数
k3 = 3.5E7; !水导轴承刚度系数
delta0 = 8.0E-3; !发电机转子不偏心时平均气隙长度
e1 = 0.2E-3 !发电机转子偏心距离
e2 = 0.3E-3; !水轮机转轮偏心距离
mu0 = 4.0*PIE-7; !空气磁导系数
kj = 5.2E0; !气隙基波磁动势系数
R = 1.2E0; !发电机转子半径
L = 0.5E0; !发电机转子长度
omega = 10.0E0; !机组旋转频率 单位是Hz
!!!!! 转轮密封力参数 !!!!!!
alpha0 = 1.0E0; !动量修正系数
rho = 1.0E0; !液体密度，单位kg/m^3
lambda = 1.0E0; ! 沿程损失系数
xi = 1.5E0; ! 局部损失系数
Rm = 2.925E0; ! 密封半径 单位m
lm = 0.43E0; ! 密封长度 单位m
v = 3.537E0; ! 液体轴向流速 单位m/s
cm = 2.5E-3; ! 密封间隙 单位m
K_1 = 25*k1 + 9*k2 + k3
K_2 = -5*k1 + 3*k2 + 3*k3
K_3 = k1 + k2 + 9*k3
B1 = e1*omega**2 * cos(2*PI*omega*T)/delta0
B2 = e1*omega**2 * sin(2*PI*omega*T)/delta0
Kxx = lambda*lm**2 *rho*PI*v**2 *Rm*(1.0+xi)/(8.0*cm**2 *(1.0+xi+lambda*lm/2.0/cm))
Kyy = Kxx
Kyx = -(Rm*omega/2.0)*(PI*rho*v*lm**2 /2.0/cm)*(1.0+xi+lambda*lm/6.0/cm)
Kxy = -Kyx ! 刚度系数
xishu = (2.0*(1+xi)+lambda*lm/4.0/cm)/(3.0*(1.0+xi)+lambda*lm/2.0/cm)
Dxx = (rho*v*lm**2 /2.0/cm)*PI*Rm*(1.0+xi+lambda*lm/6.0/cm)
Dyy = Dxx
Dyx = -(Rm*omega/2.0)*(PI*rho*alpha0*lm**3 /cm)*xishu
Dxy = -Dyx !阻尼系数
mxx = PI*Rm*rho*alpha0*lm**3 /cm*xishu;
myy = mxx ! 惯性系数
M2x = m2 + mxx
M2y = m2 + myy
D1 = (c2 + Dxx)/M2x
D2 = Dxy/M2x
D3 = (c2 + Dyy)/M2y
D4 = Dyx/M2y
Kg1 = Kxy/M2x
Kg2 = Kyx/M2y
B3 = m2*e2*omega**2 * cos(2*pi*omega*T)/(M2x*delta0)
B4 = m2*e2*omega**2 * sin(2*pi*omega*T)/(M2y*delta0)
Z1 = (1.0 - sqrt(1.0 - Y(1)**2 -Y(3)**2 )) / sqrt(Y(1)**2 + Y(3)**2)
Z2 = sqrt(Y(5)**2 + Y(7)**2) / sqrt(Y(1)**2 + Y(3)**2)
Z3 = 1.0 / Z2
!!!!! 不平衡磁拉力 !!!!!!11
!e = sqrt(Y(1)**2 + Y(3)**2)
F0 = R*L*PI*mu0*(Ij*Kj)**2 / (m1*delta0**3)
F1 = (Z1 + Z1**3 + Z1**5) / (1-Y(1)**2 - Y(3)**2)
F = F0 * F1
!!!!!!! 系统运动微分方程 !!!!!!!!!!!!!1
YPRIME(1) = Y(2)
YPRIME(2) = B1 + F*Y(1)/sqrt(Y(1)**2 + Y(3)**2) - c1/m1*Y(2) - 1.0/(16.0*m1) * (K_1+K_2*Z2) * Y(1)
YPRIME(3) = Y(4)
YPRIME(4) = B2 + F*Y(3)/sqrt(Y(1)**2 + Y(3)**2) - c1/m1*Y(4) - 1.0/(16.0*m1) * (K_1+K_2*Z2) * Y(3)
YPRIME(5) = Y(6)
YPRIME(6) = B3 - D1*Y(6) - D2*Y(8) - 1.0/(16.0*M2x) * (K_3+K_2*Z3) * Y(5) - Kg1 * Y(7)
YPRIME(7) = Y(8)
YPRIME(8) = B4 - D3*Y(8) - D2*Y(6) - 1.0/(16.0*M2y) * (K_3+K_2*Z3) * Y(7) - Kg2 * Y(5)
return
end subroutine FCN
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine FCNJ (N, T, Y, DYPDY)
integer N
real T, Y(N), DYPDY(N,*)
return
end subroutine FCNJ ! 实际上并没有用到FCNJ
[/code]
结果提示如图所示：
之前一直用matlab计算微分方程组，所以对fortran不是很了解。但Fortran的计算速度是matlab无法比拟的，所以想修改程序。麻烦大家看一下究竟是什么原因造成这种现象的（方程应该不是刚性的，否则ode45根本不可能在2.5秒之内就计算出来）。