Fortran转换matlab计算微分方程组，提示：可能为刚性问题

chunshui2003 · 发表于 2010-10-29 21:20

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

matlab的代码如下：

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);

复制代码

求解：

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

复制代码

大概计算时间为2.5秒，没有问题，非刚性。

改为Fortran程序如下：

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

复制代码

结果提示如图所示：

Fortran运行提示

之前一直用matlab计算微分方程组，所以对fortran不是很了解。但Fortran的计算速度是matlab无法比拟的，所以想修改程序。麻烦大家看一下究竟是什么原因造成这种现象的（方程应该不是刚性的，否则ode45根本不可能在2.5秒之内就计算出来）。

Jillian · 发表于 2010-10-31 05:59

本帖最后由 Jillian 于 2010-10-31 06:00 编辑

你这个imsl库是哪个版本的？好像和我用的调用不是太一样
我用的一般前面都要最一些初始化

call umach(2,nout) ！初始化计算条件

复制代码

call dset(mxparm,0.0d0,param,1) ！参数初始化
param(4)=160000
param(10)=1.0

复制代码

也没有你这里的雅克比矩阵之类的，另外函数调用的格式也不一样

call ivprk(ido,n,fcn,t,tend,tol,param,x)

复制代码

最好先找对应的帮助文件中的例子做一下看看

chunshui2003 · 发表于 2010-10-31 10:02

回复 Jillian 的帖子

Jillian，你好。我的fortran版本为CVF6.6，IMSL是随软件一起携带的4.0版（1998，应该是这样）。你所提到的初始化条件我给省略掉了，因为在实验一个简单的例子时发现是否有

call umach(2,nout)

复制代码

对计算结果几乎没有影响。经你提醒后我会再加上看一下。

param(4)=160000 !设置最大步数
param(10) = 1.0 !采用绝对误差

复制代码

这两个地方没有问题，尽管我设置了100000步也是无济于事。
而IVPRK和IVPAG这两个命令我都有仔细阅读帮助和看例子。二者区别在于：
1. IVPRK主要用来求解非刚性微分方程，尽管也能够求解刚性问题，但代价较大，计算步数过多并且有可能出现无法求解的现象
2. IVPAG可以求解刚性，也可以求解非刚性问题。
因为我实际要算的东西是刚性的（不是帖子中的例子），在matlab中已经得到了验证，用ode15s大概也要计算20秒左右，ode45根本无法求解。
所以，我先拿以前我算过的一个例子实验一下，看IVPAG这个函数是否管用。如果好使的话，我再进行下一步的操作。
结果，就出现了上述的问题。
而我刚刚接触Fortran不久，所以只能是照着葫芦画瓢，有很多细节上的问题不太清楚。毕竟，matlab和Fortran的语法结构有较大差异。如果有可能的话，也请你给予一些指点，我将非常感激。

chunshui2003 · 发表于 2010-10-31 11:09

回复 Jillian 的帖子

初始条件添加后，采用IVPRK的算法，如下：

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
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.01 ! 绝对误差
call umach(2,nout)
call sset (MXPARM, 0.0, PARAM, 1)
!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 IVPRK (IDO, N, FCN, T, TEND, TOL, PARAM, Y)
write(NOUT, '(I6, 9F12.3)') ISTEP, T, Y
end do
call IVPRK (3, N, FCN, 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 PI
save PI
real F0, F1, F, e
save F0
!real, parameter :: Ij = 500.0
real Ij
save Ij
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
PI = 3.14159
Ij = 500.0
!!!!! 运动微分方程参数 !!!!!!
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

复制代码

得到的结果：

运行结果

这个我没有看懂，希望你给些指点，是提示主程序有问题还是算法有问题呢？

Jillian · 发表于 2010-11-7 09:25

回复 4 # chunshui2003 的帖子

从提示来看，应该是

# %%%%%%%%%%% 涉及到变量的一些参数
# 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+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);

复制代码

的问题，在迭代的过程中，sqrt中出现了负号
这个问题比较头疼，你减小步长看看

风花雪月 · 发表于 2010-11-24 09:57

本帖最后由风花雪月于 2010-11-24 09:59 编辑

带有开方的这种模型在数值分析中确实很头疼，一般解决办法主要有三种途径
很容易出现DOMAIN error这类错误
1. 采取足够小的步长，不过计算代价呈几何指数增加
2. 给出足够好的初值，初值问题本身就是个难题
3. 人为判断sqrt中是否出现了负值，然后重新赋值，这种办法仅对部分问题有效

chunshui2003 · 发表于 2010-11-24 15:40

风花雪月发表于 2010-11-24 09:57

登录/注册后可看大图

带有开方的这种模型在数值分析中确实很头疼，一般解决办法主要有三种途径
很容易出现DOMAIN error这类错误 ...

谢谢风花雪月的回答，因为没人交流，这段时间又很急，所以只好暂时拿matlab算了。尽管时间长了点，但是能够运行，结果也还可以。不过，一旦闲下来还是要继续钻研fortran的，毕竟效率高很多！

chunshui2003 · 发表于 2010-11-24 16:12

感谢Rainboy，如果以后用到C或C++还要向你请教。

账号		自动登录	找回密码
密码			我要加入

[Fortran] Fortran转换matlab计算微分方程组，提示：可能为刚性问题

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

评分

评分

评分

评分

点评