哈爾濱工建大學(xué)數(shù)值大作業(yè)2014-附fortran程序

1、B班大作業(yè)要求：1. 使用統(tǒng)一封皮；2. 上交大作業(yè)內(nèi)容包含： 一 摘要 二 數(shù)學(xué)原理 三 程序設(shè)計（必須對輸入變量、輸出變量進行說明；編程無語言要求，但程序要求通過） 四 結(jié)果分析和討論 五 完成題目的體會與收獲3. 提交大作業(yè)的時間：本學(xué)期最后一次課，或考前答疑；過期不計入成績；4. 提交方式：打印版一份；或手寫大作業(yè)，但必須使用A4紙。5. 撰寫的程序需打印出來作為附錄。課 程 設(shè) 計課程名稱： 設(shè)計題目： 學(xué) 號： 姓 名： 完成時間： 題目一：非線性方程求根一 摘要非線性方程的解析解通常很難給出，因此非線性方程的數(shù)值解就尤為重要。本實驗通過使用常用的求解方法二分法和Newton法及改

2、進的Newton法處理幾個題目，分析并總結(jié)不同方法處理問題的優(yōu)缺點。觀察迭代次數(shù)，收斂速度及初值選取對迭代的影響。用Newton法計算下列方程 （1） ， 初值分別為，； （2） 其三個根分別為。當(dāng)選擇初值時給出結(jié)果并分析現(xiàn)象,當(dāng)，迭代停止。二 數(shù)學(xué)原理對于方程f(x)=0，如果f(x)是線性函數(shù)，則它的求根是很容易的。牛頓迭代法實質(zhì)上是一種線性化方法，其基本思想是將非線性方程f(x)=0逐步歸結(jié)為某種線性方程來求解。設(shè)已知方程f(x)=0有近似根xk(假定) ，將函數(shù)f(x)在點xk進行泰勒展開，有于是方程f(x)=0可近似的表示為這是個線性方程，記其根為xk+1，則xk+1的計算公式為 ，

3、k=0,1,2,這就是牛頓迭代法或簡稱牛頓法。三 程序設(shè)計（本程序由Fortran語言編制）（1）對于，按照上述數(shù)學(xué)原理，編制的程序如下 program newtonimplicit nonereal : x(0:50),fx(0:50),f1x(0:50)!分別為自變量x，函數(shù)f(x)和一階導(dǎo)數(shù)f1(x)integer : kwrite(*,*) "x(0)="read(*,*) x(0) !輸入變量：初始值x(0) open(10,file='1.txt')do k=1,50,1 fx(k)=x(k-1)*3-x(k-1)-1 f1x(k)=3*x(k-

4、1)*2-1 x(k)=x(k-1)-fx(k)/f1x(k) !牛頓法 write(*,'(I3,1x,f11.6)') k,x(k) !輸出變量：迭代次數(shù)k及x的值 write(10,'(I3,1x,f11.6)') k,x(k) if(abs(x(k)-x(k-1)<1e-6) exit !終止迭代條件end dostopend（2）對于，按照上述數(shù)學(xué)原理，編制的程序如下program newtonimplicit nonereal : x(0:50),fx(0:50),f1x(0:50)!分別為自變量x，函數(shù)f(x)和一階導(dǎo)數(shù)f1(x)intege

5、r : kwrite(*,*) "x(0)="read(*,*) x(0) !輸入變量：初始值x(0) open(10,file='1.txt')do k=1,50,1 fx(k)=x(k-1)*3+94*x(k-1)*2-389*x(k-1)+294 f1x(k)=3*x(k-1)*2+188*x(k-1)-389 x(k)=x(k-1)-fx(k)/f1x(k) !牛頓法 write(*,'(I3,1x,f11.6)') k,x(k) !輸出變量：迭代次數(shù)k及x的值 write(10,'(I3,1x,f11.6)') k,

6、x(k) if(abs(x(k)-x(k-1)<5e-6) exit !終止迭代條件end dostopend四 結(jié)果分析和討論（1）對于方程，當(dāng)初值分別為，時，所得結(jié)果如下kxk初始值1初始值2初始值30x010.450.651x11.500000-3.0121025.7915912x21.347826-2.0465173.9098533x31.325200-1.3958492.6869634x41.324718-0.9162361.9264205x51.324718-0.3545281.5097046x6-1.4622501.3501837x7-0.9701851.3253028x8

7、-0.4531211.3247189x9-2.1193701.32471810x10-1.44601211x11-0.95718212x12-0.43116813x13-1.89852714x14-1.29275915x15-0.82741716x16-0.12613717x17-1.04590918x18-0.56460119x19-14.65421020x20-9.78310721x21-6.54137022x22-4.38730123x23-2.95878924x24-2.01102125x25-1.37128326x26-0.89570027x27-0.31076928x28-1.32

8、340729x29-0.85459830x30-0.20847031x31-1.12909032x32-0.66518233x331.25644434x341.32950635x351.32473936x361.32471837x371.324718結(jié)果分析與討論：從計算結(jié)果可以看到，當(dāng)取的初始值不同時，雖然均得到了近似解x*=1.324718，但收斂速度明顯不同。當(dāng)初始值x0充分接近于方程的單根時，可保證迭代序列快速收斂。在本例中，初始值1、0.65和0.45距方程的單根越來越遠，故收斂速度越來越慢。（2）對于方程，當(dāng)初始值x0=2時計算結(jié)果如下kxk初始值0x021x1-98.000000

9、2x2-98.000000結(jié)果分析與討論：牛頓法有明顯的幾何解釋，方程f(x)=0的根x*可解釋為曲線y=f(x)與x軸的交點的橫坐標。設(shè)xk是根x*的某個近似值，過曲線y=f(x)上橫坐標為xk的點Pk引曲線y=f(x)的切線，并將該切線與x軸的交點坐標xk+1作為x*的新的近似值。在本例中，當(dāng)初始值x0=2時，在這個點的切線方程與x軸的交點恰為方程的一個根x=-98,故迭代了兩次就得到了結(jié)果。五 完成題目的體會與收獲 對于牛頓法，當(dāng)初始值x0充分接近于方程的單根時，可保證迭代序列快速收斂。當(dāng)方程有不止一個根時，所得結(jié)果與初始值的選取有關(guān) 。 題目二：線性方程組求解一 摘要對于實際的工程問題

10、，很多問題歸結(jié)為線性方程組的求解。本實驗通過實際題目掌握求解線性方程組的數(shù)值解法，直接法或間接法。有一平面機構(gòu)如圖所示，該機構(gòu)共有13條梁（圖中標號的線段）由8個鉸接點（圖中標號的圈）聯(lián)結(jié)在一起。上述結(jié)構(gòu)的1號鉸接點完全固定，8號鉸接點豎立方向固定，并在2號、5號和6號鉸接點，分別有如圖所示的10噸、15噸和20噸的負載，在靜平衡的條件下，任何一個鉸接點上水平和豎立方向受力都是平衡的，以此計算每個梁的受力情況。7865434813579111221261013101520 令，假設(shè)為各個梁上的受力，例如對8號鉸接點有對5號鉸接點，則有針對各個鉸接點，列出方程并求出各個梁上的受力。二 數(shù)學(xué)原理高

11、斯列主元消去法：方法說明（以4階為例）：第1步消元在增廣矩陣（A，b）第一列中找到絕對值最大的元素，將其所在行與第一行交換，再對（A，b）做初等行變換使原方程組轉(zhuǎn)化為如下形式：第2步消元在增廣矩陣（A，b）中的第二列中（從第二行開始）找到絕對值最大的元素，將其所在行與第二行交換，再對（A，b）做初等行變換使原方程組轉(zhuǎn)化為：第3步消元在增廣矩陣（A，b）中的第三列中（從第三行開始）找到絕對值最大的元素，將其所在行與第二行交換，再對（A，b）做初等行變換使原方程組轉(zhuǎn)化為：按x4 à x3à x2à x1 的順序回代求解出方程組的解。針對各個鉸接點列方程：對2號鉸接點有

12、對3號鉸接點有對4號鉸接點有對5號鉸接點有對6號鉸接點有對7號鉸接點有對8號鉸接點有三 程序設(shè)計（本程序由Fortran語言編制）program mainimplicit noneinteger,parameter : n=13 ！輸入量：系數(shù)陣的維數(shù)real : js(n)dimension : a(n,n),b(n),x(n)double precision a,b,xreal,parameter : m=0.7071 !m代表= integer : i,jlogical ldata(a(i,j),j=1,13),i=1,13) / m,0.,1.,0.,m,0.,0.,0.,0.,0.,

13、0.,0.,0., & ！輸入量：方程的系數(shù)陣 0.,1.,0.,0.,0.,-1.,0.,0.,0.,0.,0.,0.,0., & 0.,0.,1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., & 0.,0.,0.,1.,0.,0.,0.,-1.,0.,0.,0.,0.,0., & m,0.,0.,-1.,-m,0.,0.,0.,0.,0.,0.,0.,0., & 0.,0.,0.,0.,m,1.,0.,0.,-m,-1.,0.,0.,0.,& 0.,0.,0.,0.,0.,0.,1.,0.,0.,0.,0.,0.,0.,&a

14、mp;0.,0.,0.,0.,0.,0.,0.,1.,m,0.,0.,-m,0.,& 0.,0.,0.,0.,m,0.,1.,0.,m,0.,0.,0.,0.,& 0.,0.,0.,0.,0.,0.,0.,0.,0.,1.,0.,0.,-1.,& 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,1.,0.,0.,& 0.,0.,0.,0.,0.,0.,0.,0.,m,0.,1.,m,0.,& 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,m,1. / b=0. ！輸入量：方程的右邊項 b(3)=10. b(9)=15. b

15、(11)=20.write(*,"(13(13(f5.2,1x),/)") (a(i,j),j=1,13),i=1,13) ！輸出矩陣call agaus(a,b,n,x,l,js) if (l/=0) then write(*,"(3(5f8.2,/)") x(:) ！輸出結(jié)果end ifstopendsubroutine agaus(a,b,n,x,l,js) dimension a(n,n),x(n),b(n),js(n) double precision a,b,x,t l=1 !邏輯變量 do k=1,n-1 d=0.0 do i=k,n do

16、 j=k,n if (abs(a(i,j)>d) then d=abs(a(i,j) js(k)=j is=i end if end do end do !把行絕對值最大的元素換到主元位置 if (d+1.0=1.0) then l=0 else !最大元素為0無解 if(js(k)/=k) then do i=1,n t=a(i,k) a(i,k)=a(i,js(k) a(i,js(k)=t end do !最大元素不在K行，K行 end if if(is/=k) then do j=k,n t=a(k,j) a(k,j)=a(is,j) a(is,j)=t !交換到K列 end do

17、 t=b(k) b(k)=b(is) b(is)=t end if !最大元素在主對角線上 end if !消去 if (l=0) then write(*,100) return end if do j=k+1,n a(k,j)=a(k,j)/a(k,k) end do b(k)=b(k)/a(k,k) !求三角矩陣 do i=k+1,n do j=k+1,n a(i,j)=a(i,j)-a(i,k)*a(k,j) end do b(i)=b(i)-a(i,k)*b(k) end do end do if (abs(a(n,n)+1.0=1.0) then l=0 write(*,100)

18、return end if x(n)=b(n)/a(n,n) do i=n-1,1,-1 t=0.0 do j=i+1,n t=t+a(i,j)*x(j) end do x(i)=b(i)-t end do100 format(1x,'fail') js(n)=n do k=n,1,-1 if (js(k)/=k) then t=x(k) x(k)=x(js(k) x(js(k)=t end if end do return end四 結(jié)果分析和討論由程序計算的各個桿的受力如下：f1f2f3f4f5f6f7-28.2820.0010.00-30.0014.1420.000.00

19、f8f9f10f11f12f13-30.007.0725.0020.00-35.3625.00結(jié)果分析與討論：列方程時假定各桿均受拉力，得到的結(jié)果有負值，說明力與假設(shè)方向相反。五 完成題目的體會與收獲高斯消去法雖然能夠執(zhí)行，但隨便用akk(k)作為主元，有時會擴大誤差，導(dǎo)致結(jié)果不可靠，甚至嚴重失真，而高斯列主元消去法則不會有這種情況發(fā)生。最初采用的是高斯-賽德爾迭代法解此矩陣，但是發(fā)現(xiàn)很難得到收斂解。因為必須滿足迭代條件才可以，而找到滿足條件的迭代矩陣非常困難，故最終選取了沒有收斂限制的直接法解此矩陣。附錄題目一程序：（1）program newtonimplicit nonereal : x

20、(0:50),fx(0:50),f1x(0:50)!分別為自變量x，函數(shù)f(x)和一階導(dǎo)數(shù)f1(x)integer : kwrite(*,*) "x(0)="read(*,*) x(0) !輸入變量：初始值x(0) open(10,file='1.txt')do k=1,50,1 fx(k)=x(k-1)*3-x(k-1)-1 f1x(k)=3*x(k-1)*2-1 x(k)=x(k-1)-fx(k)/f1x(k) !牛頓法 write(*,'(I3,1x,f11.6)') k,x(k) !輸出變量：迭代次數(shù)k及x的值 write(10,&#

21、39;(I3,1x,f11.6)') k,x(k) if(abs(x(k)-x(k-1)<1e-6) exit !終止迭代條件end dostopend（2）program newtonimplicit nonereal : x(0:50),fx(0:50),f1x(0:50)!分別為自變量x，函數(shù)f(x)和一階導(dǎo)數(shù)f1(x)integer : kwrite(*,*) "x(0)="read(*,*) x(0) !輸入變量：初始值x(0) open(10,file='1.txt')do k=1,50,1 fx(k)=x(k-1)*3+94*x(

22、k-1)*2-389*x(k-1)+294 f1x(k)=3*x(k-1)*2+188*x(k-1)-389 x(k)=x(k-1)-fx(k)/f1x(k) !牛頓法 write(*,'(I3,1x,f11.6)') k,x(k) !輸出變量：迭代次數(shù)k及x的值 write(10,'(I3,1x,f11.6)') k,x(k) if(abs(x(k)-x(k-1)<5e-6) exit !終止迭代條件end dostopend題目二程序：program mainimplicit noneinteger,parameter : n=13 ！輸入量：系數(shù)陣的

23、維數(shù)real : js(n)dimension : a(n,n),b(n),x(n)double precision a,b,xreal,parameter : m=0.7071 !m代表= integer : i,jlogical ldata(a(i,j),j=1,13),i=1,13) / m,0.,1.,0.,m,0.,0.,0.,0.,0.,0.,0.,0., & ！輸入量：方程的系數(shù)陣 0.,1.,0.,0.,0.,-1.,0.,0.,0.,0.,0.,0.,0., & 0.,0.,1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., & 0.,

24、0.,0.,1.,0.,0.,0.,-1.,0.,0.,0.,0.,0., & m,0.,0.,-1.,-m,0.,0.,0.,0.,0.,0.,0.,0., & 0.,0.,0.,0.,m,1.,0.,0.,-m,-1.,0.,0.,0.,& 0.,0.,0.,0.,0.,0.,1.,0.,0.,0.,0.,0.,0.,&0.,0.,0.,0.,0.,0.,0.,1.,m,0.,0.,-m,0.,& 0.,0.,0.,0.,m,0.,1.,0.,m,0.,0.,0.,0.,& 0.,0.,0.,0.,0.,0.,0.,0.,0.,1.,0.,0

25、.,-1.,& 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,1.,0.,0.,& 0.,0.,0.,0.,0.,0.,0.,0.,m,0.,1.,m,0.,& 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,m,1. / b=0. ！輸入量：方程的右邊項 b(3)=10. b(9)=15. b(11)=20.write(*,"(13(13(f5.2,1x),/)") (a(i,j),j=1,13),i=1,13) ！輸出矩陣call agaus(a,b,n,x,l,js) if (l/=0) then write(*,"(3(5f8.2,/)") x(:) ！輸出結(jié)果end ifstopendsubroutine agaus(a,b,n,x,l,js) dimension a(n,n),x(n),b(n),js(n) double precision a,b,x,t l=1 !邏輯變量 do k