下载文档原格式
数值计算基础目录实验一直接法解线性方程组的 (2)实验二插值方法 (11)实验三数值积分 (5)实验四常微分方程的数值解 (7)实验五迭代法解线性方程组与非线性方程 (9)实验一 直接法解线性方程组一、实验目的掌握全选主元消去法与高斯-塞德尔法解线性方程组。
二、实验内容分别写出Guass 列选主元消去法与追赶法的算法，编写程序上机调试出结果，要求所编程序适用于任何一解线性方程组问题，即能解决这一类问题，而不是某一个问题。
实验中以下列数据验证程序的正确性。
1、用Guass 列选主元消去法求解方程组⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--5.58.37.33.47.11.85.16.93.51.53.25.2321x x x2、用追赶法求解方程组⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-----000010210000210000210000210000254321x x x x x 三、实验仪器设备与材料主流微型计算机四、实验原理1、Guass 列选主元消去法 对于AX =B1）、消元过程：将（A|B ）进行变换为)~|~(B A ，其中A ~是上三角矩阵。
即：⎪⎪⎪⎪⎪⎭⎫⎝⎛→⎪⎪⎪⎪⎪⎭⎫ ⎝⎛n nn nnn nnn n nn b a b a b a a b a a a b a a a b a a a0010122111221222221111211 k 从1到n-1a 、 列选主元选取第k 列中绝对值最大元素ik ni k a ≤≤max 作为主元。
b 、 换行ik ij kj b b n k j a a ⇔+=⇔,,1,c 、 归一化kkk k kj kk kj b a b n k j a a a ⇒+=⇒/,,1,/d 、 消元nk i b b a b n k j n k i a a a a i k ik i ij kj ik ij ,,1,,,1;,,1, +=⇒-+=+=⇒-2）、回代过程：由)~|~(B A 解出11,,,x x x n n -。
(1) 采用fortran 进行计算的好处简单的计算，例如2/1-+=d k i s σσ，可以直接在origin 等里面计算；复杂的计算，例如 2220)()()(geo the or s σσσσσ∆+∆+∆+= 其中 Lbor μϕσ=∆，ρμησb the =∆，D b f m geo /εμβσ=∆又其中 2/14⎪⎪⎭⎫⎝⎛=fD L π，)1(12f D b fT -∆∆=αρ，怎么算？？(2) fortran 语言 (以以上计算为例)REAL EM,EP,YSM,F,VM,YS (定义参量 REAL ：实数，赋值时需要加.) REAL ROU REAL D,B,UM,DETAALFA,DETA T,YSINCROR,YSINCRTHE,YSINCRGEO,YSMINC This project is used to calculate the YIELD STRENGTH of AL-AMORPHOUSPARTICLES SL METHOD (说明这个程序的目的 开头C 决定了这一行是不运行的)OPEN(4,FILE='YS.DA T', STATUS='UNKNOWN') (建立数据文件.DAT ，存储数据) OPEN(5,FILE='YSINCROR.DA T', STATUS='UNKNOWN') OPEN(6,FILE='YSINCRTHE.DA T', STATUS='UNKNOWN')EM=50. (基体M 弹性模量赋值) VM=0.33 (基体M 泊松比赋值) UM=28. (基体M 剪切模量赋值) EP=95. (第二相P 弹性模量赋值)B=0.286 (基体Burges 矢量赋值) (注：为减小计算量，以上赋值均省去单位，最后一起核算)YSM=113. (基体屈服强度0σ赋值)D=100. (第二相颗粒直径赋值) DETAALFA=1.2E-5 (热膨胀系数差α∆赋值) DETAT=345. (热膨胀系数差T ∆赋值)DO 30 F=0.0001,0.8,0.0999 (DO CONTINUE 循环语句，对颗粒含量F 连续赋值) L=0.5*D*(4.188/F)**0.334 (求第二相颗粒间距L ) YSINCROR=UM*B/L (计算or σ∆)ROU=12.*DETAALFA*DETAT*F/(B*D*(1.-F)) (计算位错密度ρ) YSINCRTHE=1.4*UM*B*SQRT(1.E3*ROU) (计算the σ∆)YSINCRGEO=UM*SQRT(1.E-3*F*0.002*B/D) (计算geo σ∆) YS=YSM+SQRT(YSINCROR**2.+YSINCRTHE**2.+YSINCRGEO**2.) (计算s σ) WRITE (4,*)100.*F,YS (输出 F, s σ结算结果 )WRITE (5,*)100.*F,YSINCROR (输出 F, or σ∆结算结果 )WRITE (6,*)100.*F,YSINCRGEO,YSINCRTHE (输出 F, geo σ∆, the σ∆结算结果 )30 CONTINUE END (结束)(2) fortran 计算界面(i) compile 编辑，查错(ii) built建立文件关联(iii) go 运行计算程序存盘时必须是.f查看计算结果：(3) fortran 计算结果的输出 (origin)100120140160180Y i e l d s t r e n g t h (M P a )Volume fraction of particles (%)D=20 μm D=50 μm D=100 μm点击openAll filesF σS点击后找到相应文件附1：Al-Cu合金系析出强化计算程序涉及到的理论及公式请参考Liu et al., Materials Science and Engineering A 344(2003)113REAL R0,D,B,C0,CE,CP,P0,PM,W0,WM,G,B1,TM,FM,A0,FV1,FV2,FV,K,K1REAL D1,G1,BTIME,BTM,R0D1,T1,T2,F,G11,G12,G21,G22,A1,A2,R1,R,NINTEGER T,NOMBER,TIMEOPEN(4,FILE='AL-CU52-SHIYAN.dat')OPEN(3,FILE='AL-CU52-SHIYAN-N.OUT')R=8.31K=27.0READ(4,*) TM,A0,T,W0,NOMBERWRITE(3,*)NOMBERWRITE(3,*)KP0=8.93PM=2.6989CP=0.33WM=1.-W0C0=W0*PM/(W0*PM+WM*P0)CE=13.176*EXP(-41680./(R*T))FM=C0-CEB=0.625*FM/(K*(0.33-CE))FV1=CP*LOG(C0/CE)FV2=(1.-CP)*LOG((1.-C0)/(1.-CE))FV=R*T*(FV1+FV2)R0=8.7E3/FVD=6.0E-6*EXP(-119691./(R*T))D1=D*1.0E20BTM=(B*TM)**1.5R0D1=(R0*D1)**2.0K1=K**2.0G11=9.12E7*K1*BTM*SQRT(D)G12=R0D1*C0*TMG1=3.0*FM/(G11*G12)G=-R*T*LOG(G1)write(3,*)G/1000.DO 10 TIME=100,800,10T2=TIME/100.0T1=10.0**(T2)R1=0.666667*K*SQRT(D*B*T1)R=1.0E9*R1BTIME=(B*T1)**1.5G21=9.12E7*K1*BTIME*SQRT(D)G22=R0D1*C0*T1F=G21*G22*G1N=1.0E10*F*K/(2.5*6.28*R**3.)B1=0.286IF (0.079*R.LE.B1) THENA2=A0WRITE(3,300) T1,NELSE IF (T1.LE.TM) THENA1=1.61355E3*SQRT(K)/RA2=A0+A1*(SQRT(F)+0.75*SQRT(K)*F+0.14*K*(F)**1.5)*LOG(0.079*R/B1)WRITE(3,300) T1,NELSEF=3.0*FMA1=1.61355E3*SQRT(K)/RA2=A0+A1*(SQRT(F)+0.75*SQRT(K)*F+0.14*K*(F)**1.5)*LOG(0.079*R/B1)N=1.0E10*F*K/(2.5*6.28*R**3.)WRITE(3,300) T1,N300 FORMAT(1X,F18.6,2X,F18.6)END IF10 CONTINUECLOSE(3)CLOSE(4)END得到的结果举例：附2：复合材料应力应变曲线计算涉及到的理论及公式请参考Liu et al., Acta Materialia 57(2009)2029REAL KM,KP,KX,UMS,UM,UP,UX,VVM,VVPREAL STRESSM,STRAINM,STRAINMM,STRESSX,STRAINX,STRESSM0REAL STRESSP,STRESSPCRI,DIAM,FCRI,M,K,VOLFRREAL YS,EMS,EM,EP,EX,ALFA,N,NI,F,FF,FM,F1REAL A1,A2,A3,B1,B2,C1,C2,D1,D2,YIELDREAL STRESSX1R,STRAINXR,DETAINTEGER IC This project is used to calculate the strain-stress curves of compositesOPEN(4,FILE='STRAIN-STRESS-60%-70.DAT', STATUS='UNKNOWN') OPEN(5,FILE='STRAIN-STRESS-MATRIX.DAT', STA TUS='UNKNOWN') OPEN(6,FILE='YIELD.DA T', STATUS='UNKNOWN')OPEN(7,FILE='VOLFRA.DAT', STATUS='UNKNOWN')ALFA=0.42857N=0.03NI=1./NF=0.10F1=FFM=1.-FFF=F/FMDIAM=25.K=1200.C K=5600.STRESSPCRI=K/SQRT(DIAM)M=5.8EM=70.EP=115.EX=EM*(1.-F)+EP*FVVM=0.33VVP=0.22KM=EMUM=0.5*KM/(1.+VVM)KP=EPUP=0.5*KP/(1.+VVP)STRESSM0=205.STRESSM=0.VOLFR1=0.STRAINXR=0.STRESSXR=0.WRITE(4,*) 0,0WRITE(5,*) 0,0DO 30 I=1,400STRESSM=STRESSM+1.A1=STRESSM/STRESSM0A2=N/(1.-N)A3=1.+ALFA*(A1)**A2STRAINM=STRESSM/EM+ALFA*(STRESSM0/EM)*(A1)**NIEMS=EM/A3UMS=EMS/(3.-EMS/(3.*KM))B1=UMS*(9.*KM+8.*UMS)B2=6.*(KM+2.*UMS)YS=B1/B2YX=YSC1=3.*KP*KM*(1.+FF)+4.*UX*(KM+KP)C2=4.*UX*(1.+FF)+3.*(KP+FF*KM)KX=C1/C2D1=UMS*UP*(1.+FF)+YX*(UMS+FF*UP)D2=YX*(1.+FF)+UP+FF*UMSUX=D1/D2STRAINX=STRAINM*(3.*KM+4.*UMS)/(3.*KX+4.*UMS) STRESSX1=STRESSM*(KX/KM)*(3.*KM+4.*UMS)/(3.*KX+4.*UMS) STRESSX2=1.E3*EX*(1.E-3*STRAINX-0.002)STRESSX=STRESSX1-STRESSX2STRESSP=STRESSM*(KP/KM)*(3.*KM+4.*UMS)/(3.*KP+4.*UMS) VOLFR=1.-EXP(-(STRESSP/STRESSPCRI)**M)F=F1-F1*VOLFRFF=F/FMEM=70.*FM/(FM+F1*VOLFR)N=0.03*FM/(FM+F1*VOLFR)N1=1./NDETA=1.E3*(STRESSX1-STRESSXR)/(STRAINX-STRAINXR) WRITE (4,*)1.E-1*STRAINX,STRESSX1,DETA WRITE (5,*)1.E-1*STRAINM,STRESSM,STRESSX1,STRESSXR WRITE (7,*)1.E-1*STRAINX,VOLFR*100. STRAINXR=STRAINX STRESSXR=STRESSX1 IF(ABS(STRESSX).LT.5.) THEN YIELD=STRESSX1 ELSE ENDIF30 CONTINUEWRITE (6,*)YIELD END得到的结果举例：St r e s s (M P a )Strain (%)S t r e s s (M P a )Strain (%)。
