哈尔滨工程大学数值分析大作业附fortran程序.docx

1、哈尔滨工程大学数值分析大作业附fortran程序B班大作业要求：1. 使用统一封皮；2. 上交大作业内容包含： 一 摘要 二 数学原理 三 程序设计（必须对输入变量、输出变量进行说明；编程无语言要求，但程序要求通过） 四 结果分析和讨论 五 完成题目的体会与收获3. 提交大作业的时间：本学期最后一次课，或考前答疑；过期不计入成绩；4. 提交方式：打印版一份；或手写大作业，但必须使用A4纸。5. 撰写的程序需打印出来作为附录。课 程 设 计课程名称： 设计题目： 学 号： 姓 名： 完成时间： 题目一：非线性方程求根一 摘要非线性方程的解析解通常很难给出，因此非线性方程的数值解就尤为重要。本实验

2、通过使用常用的求解方法二分法和Newton法及改进的Newton法处理几个题目，分析并总结不同方法处理问题的优缺点。观察迭代次数，收敛速度及初值选取对迭代的影响。用Newton法计算下列方程 （1） ， 初值分别为，； （2） 其三个根分别为。当选择初值时给出结果并分析现象,当，迭代停止。二 数学原理对于方程f(x)=0，如果f(x)是线性函数，则它的求根是很容易的。牛顿迭代法实质上是一种线性化方法，其基本思想是将非线性方程f(x)=0逐步归结为某种线性方程来求解。设已知方程f(x)=0有近似根xk(假定) ，将函数f(x)在点xk进行泰勒展开，有于是方程f(x)=0可近似的表示为这是个线性方

3、程，记其根为xk+1，则xk+1的计算公式为 ，k=0,1,2,这就是牛顿迭代法或简称牛顿法。三 程序设计（本程序由Fortran语言编制）（1）对于，按照上述数学原理，编制的程序如下 program newtonimplicit nonereal : x(0:50),fx(0:50),f1x(0:50)!分别为自变量x，函数f(x)和一阶导数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(

4、k-1)*2-1 x(k)=x(k-1)-fx(k)/f1x(k) !牛顿法 write(*,(I3,1x,f11.6) k,x(k) !输出变量：迭代次数k及x的值 write(10,(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，函数f(x)和一阶导数f1(x)integer : kwrite(*,*) x(0)=rea

5、d(*,*) 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) !输出变量：迭代次数k及x的值 write(10,(I3,1x,f11.6) k,x(k) if(abs(x(k)-x(k-1)d) then d=abs(a(i,j) js(k)=j is=i end if end

6、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 t=b(k) b(k)=b(is) b(is)=t end if !最大元素在主对角线上 end if !消去 if (l=0) then

7、 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) 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,

8、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四 结果分析和讨论由程序计算的各个杆的受力如下：f1f2f3f4f5f6f7-28.2820.0010.00-30.0014.1420.000.00f8f9f10f11f12f13-30.007.0725.0020.00-35.3625.00结果分析与讨论：列方程时假定各杆均受拉力，得到的结果有负值，说明力

9、与假设方向相反。五 完成题目的体会与收获高斯消去法虽然能够执行，但随便用akk(k)作为主元，有时会扩大误差，导致结果不可靠，甚至严重失真，而高斯列主元消去法则不会有这种情况发生。最初采用的是高斯-赛德尔迭代法解此矩阵，但是发现很难得到收敛解。因为必须满足迭代条件才可以，而找到满足条件的迭代矩阵非常困难，故最终选取了没有收敛限制的直接法解此矩阵。附录题目一程序：（1）program newtonimplicit nonereal : x(0:50),fx(0:50),f1x(0:50)!分别为自变量x，函数f(x)和一阶导数f1(x)integer : kwrite(*,*) x(0)=rea

10、d(*,*) 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) !输出变量：迭代次数k及x的值 write(10,(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),

11、fx(0:50),f1x(0:50)!分别为自变量x，函数f(x)和一阶导数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(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) !输出变量：迭代次数k及x的值 write(10,(I3,1x,f11.

12、6) k,x(k) if(abs(x(k)-x(k-1)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

13、 do 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) 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