关于fortran编程问题

caizi2008

最近用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
complex  EM,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

心灯

问题太笼统，程序求解什么方程，出现什么错误，好多程序背后的信息，都没有说明，别人没法下手。

caizi2008

就是求解f(n)的表达式,其中xa(n)是未知数
大侠运行下吧 其实具体出现的什么错误我也不大清楚
呵呵 菜鸟
问题主要出现在 call......那行

suffer

应该是调用有问题，自己参考例子修改一下吧

caizi2008

非常感谢