关于fortran编程问题
最近用fortran程序解方程,由于刚刚用这个软件,好多地方都不懂呢我用这个程序解方程,可是出现错误,请大家给予指点,说明下:我要得到的9个未知数结果可能都是复数
原程序代码:
program main
use imsl
implicit none
external FCN
complex, parameter::errrel=cmplx(0.0001,0.0001)
integer,parameter::n=9!!!!!!!!!!!!!!
integer,parameter::itmax=100
complex::xguess(n)=(/cmplx(0.0,0.0),cmplx(0.0,0.0),cmplx(0.0,0.0),cmplx(0.0,0.0),cmplx(0.0,0.0),cmplx(0.0,0.0),cmplx(0.0,0.0),cmplx(0.0,0.0),cmplx(0.0,0.0)/) !(/0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0/)
complex x(n),fnorm
call neqnf(fcn,errrel,n,itmax,xguess,x,fnorm)
write(*,*)x
end
!!!!!!!!!!!!!!
subroutine fcn(xa,f,n)
implicit none
integer n
complex,target::xa(n)
complex f(n)
real,pointer::x,y,z
real xigmaM,xigmaI,xigmac,ra,dertaa,fI
complexEM,EI,Ec,lamdaM,lamdaI,lamdac,KM,KI,Kc,miuM,miuI,miuc,&
&alphaM,alphaI,alphac,beitaM,beitaI,beitac,faic,fsaic,&
&J2,J3,J4,J5,J6,J7,J8,J9,J10,K1,K2,K3,K4,K5,K6
ra=0.5 !%%% 粒子空腔填充物半径,已知
dertaa=0.001!%%% 粒子包覆层厚度,已知
fI=0.005 !%%% 散射粒子孔隙率,已知
EM=cmplx(6e8,0.45*6e8) !%%% 基材弹性模量,已知
EI=cmplx(5.4e8,0.32*5.4e8) !%%% 散射粒子中空材料弹性模量,已知
Ec=cmplx(5e8,0.7*5e8) !%%% 包覆层材料弹性模量,已知
xigmaM=0.45 !%%% 基材泊松比,已知
xigmaI=0.32 !%%% 散射粒子中空材料泊松比,已知
xigmac=0.37 !%%% 包覆层材料泊松比,已知
lamdaM=EM*xigmaM/((1.+xigmaM)*(1.-2.*xigmaM))!%%%%%%%基材拉密常数,可算
lamdaI=EI*xigmaI/((1.+xigmaI)*(1.-2.*xigmaI))!%%%%%%%散射粒子中空材料拉密常数,可算
lamdac=Ec*xigmac/((1.+xigmac)*(1.-2.*xigmac))!%%%%%%%包覆层材料拉密常数,可算
miuM=EM/(2.*(1.+xigmaM)) !%%%%%%%基材切变模量,可算
miuI=EI/(2.*(1.+xigmaI)) !%%%%%%%散射粒子中空材料模量,可算
miuc=Ec/(2.*(1.+xigmac)) !%%%%%%%包覆层材料切变模量,可算
KM=EM/(3.*(1.-2.*xigmaM)) !%%%%%%%基材体积压缩模量,可算
KI=EI/(3.*(1.-2.*xigmaI)) !%%%%%%%散射粒子中空材料体积压缩模量,可算
Kc=Ec/(3.*(1.-2.*xigmac)) !%%%%%%%包覆层材料体积压缩模量,可算
alphaM=1/3*(1.+xigmaM)/(1.-xigmaM) !%%%%%,可算
alphac=1/3*(1.+xigmac)/(1.-xigmac) !%%%%%,可算
alphaI=1/3*(1.+xigmaI)/(1.-xigmaI) !%%%%%,可算
beitaM=2/15*(4.-5.*xigmaM)/(1.-xigmaM)!%%%%%,可算
beitaI=2/15*(4.-5.*xigmaI)/(1.-xigmaI)!%%%%%,可算
beitac=2/15*(4.-5.*xigmac)/(1.-xigmac)!%%%%%,可算
faic=(Kc+alphac*(KI-Kc))/Kc !%%%%%,可算
fsaic=(miuc+beitac*(miuI-miuc))/miuc !%%%%%,可算
!!!!!!!!!!!!!!!!!!!!!!!!J\K是为了简化方程而得到的中间参数
J2=3.*dertaa/ra*(miuI-miuc)/miuc
J3=3.*dertaa/ra*(KI-Kc)/Kc
J4=Kc*(1.-alphac)
J5=miuc*(1.-beitac)
J6=fI*(lamdaI-lamdaM)
J7=fI*2/3*(miuI-miuM)
J8=fI*3.*dertaa/ra*(lamdac-lamdaM)*faic
J9=fI*3.*dertaa/ra*2/3*(miuc-miuM)
J10=fI*((miuI-miuM)+3.*dertaa/ra*(miuc-miuM)*fsaic)
K1=J6+J7+J8+J9
K2=J7+J9*fsaic
K3=1.+J3*alphac
K4=KI-J3*J4
K5=1.+J2*beitac
K6=miuI-J2*J5
!!!miue-xa(1);Ee-xa(2);xigmae-xa(3);Ke-xa(4);lamdae-xa(5);alphae-xa(6);beitae-xa(7);A-xa(8);B-xa(9)
f(1)=xa(1)-xa(2)/(2*(1+xa(3)))
f(2)=xa(4)-xa(2)/(3*(1-2*xa(3)))
f(3)=xa(5)-xa(2)*xa(3)/((1+xa(3))*(1-2*xa(3)))
f(4)=xa(6)-1/3*(1+xa(3))/(1-xa(3))
f(5)=xa(7)-2/15*(4-5*xa(3))/(1-xa(3))
f(6)=xa(8)-xa(4)/(xa(4)*(1-xa(6))*K3+xa(6)*K4)
f(7)=xa(9)-xa(1)/(xa(1)*(1-xa(7))*K5+xa(7)*K6)
f(8)=xa(5)-lamdaM+K1*xa(8)-K2*xa(9)
f(9)=xa(1)-miuM+xa(9)*J10
!f(1)=xa(1)*xa(1)+xa(2)*xa(2)+xa(3)*xa(3)-3 !!!!!!!!!!
!f(2)=xa(1)*xa(2)+xa(2)*xa(3)+xa(1)*xa(3)-3 !!!!!!!!!!
!f(3)=exp(xa(1))+exp(xa(2))+exp(xa(3))-3*exp(1.)!!!!!!!!!!
return
end subroutine
回复 #1 caizi2008 的帖子
问题太笼统，程序求解什么方程，出现什么错误，好多程序背后的信息，都没有说明，别人没法下手。 就是求解f(n)的表达式,其中xa(n)是未知数大侠运行下吧 其实具体出现的什么错误我也不大清楚
呵呵 菜鸟
问题主要出现在 call......那行 应该是调用有问题，自己参考例子修改一下吧
非常感谢
页:
[1]