兰州交通大学fortran课程设计之欧阳育创编.docx
《兰州交通大学fortran课程设计之欧阳育创编.docx》由会员分享,可在线阅读,更多相关《兰州交通大学fortran课程设计之欧阳育创编.docx(17页珍藏版)》请在冰豆网上搜索。
兰州交通大学fortran课程设计之欧阳育创编
时间:
2021.02.04
创作:
欧阳育
1求一元方程的根和求定积分
1.1求一元方程的根
作业要求
1、采用函数子程序定义一元方程;
2、程序选择以下三种方法求该方程的根;
METHOD=1牛顿迭代法
METHOD=2二分法
1METHOD=3弦截法
3、对于不同的近似算法分别编写子程序,精度要求10-6。
本题用二分法、弦解法和牛顿迭代法求x3-2x2+7x+4=0的根来编写程序求解。
二分法
二分法基本思路:
现任取两个值x1和x2,使得f(x1)*f(x2)<0,也就是f(x1)和f(x2)必须异号。
这才能保证在[x1,x2]区间有解,即存在一个x使得f(x)=0。
令x=(x1+x2)/2,如果f(x)=0,就找到了这个解,计算完成。
由于f(x)是一个实型数据,所以在判断f(x)是否等于0时,是通过判断|f(x)|是否小于一个很小的数ε,如果是就认为f(x)=0。
若f(x)不等于0,判断如果f(x1)和f(x)异号,就说明解在[x1,x]区间,就以x1,x为新的取值重复步骤
(2),这时用x代替否则x2,否则反之,直到找到满足条件的解为止。
程序编写如下:
programlt12_1
realx1,x2,x
realbisect,func!
对要调用的子程序作说明
do!
输入x1和x2直到f(x1)和f(x2)异号为止
print*,'输入x1,x2的值:
'
read*,x1,x2
if(func(x1)*func(x2)<0.0)exit
print*,'不正确的输入!
'
enddo
x=bisect(x1,x2)!
调用二分法求解函数
print10,'x=',x!
输出计算结果
10format(a,f15.7)
end
realfunctionbisect(x1,x2)!
二分法求解函数
realx1,x2,x,f1,f2,fx
x=(x1+x2)/2.0
fx=func(x)
dowhile(abs(fx)>1e-6)
f1=func(x1)
if(f1*fx<0)then
x2=x
else
x1=x
endif
x=(x1+x2)/2.0
fx=func(x)
enddo
bisect=x
end
functionfunc(x)!
需要求解的函数
realx
func=x**3-2*x**2+7*x+4
end
运行结果:
1.1.1弦截法
弦截法的基本思路:
现任取两个值x1和x2,使得f(x1)*f(x2)<0。
(1)做一条通过(x1,f(x1))和(x2,f(x2))两点的直线,这条直线与x轴的交点为x。
可用以下公式求出
X=x2-(x2-x1)*f(x2)/(f(x1)-f(x2)),
(2)代入函数求得f(x),判断|f(x)|是否小于一个很小的数ε,如果是就认为f(x)=0。
(3)否则,判断如果f(x1)和f(x)异号,就说明解在[x1,x]区间,就以x1,x为新的取值重复步骤
(2),否则反之,然后以同样的办法再进一步缩小范围,直到|f(x)|<ε。
程序编写如下:
realx1,x2,x
realsecant,func!
对要调用的子程序作说明
do!
输入x1和x2直到f(x1)和f(x2)异号为止
print*,'输入x1,x2的值'
read*,x1,x2
if(func(x1)*func(x2)<0)exit
print*,'不正确的取值'
enddo
x=secant(x1,x2)!
调用弦截法求解函数
print10,'x=',x!
输出计算结果
10format(a,f15.7)
End
realfunctionsecant(x1,x2)!
弦截法求解函数
implicitnone
realx1,x2,x,f1,f2,fx
realfunc
x=x2-(x2-x1)/(func(x2)-func(x1))*func(x2)
fx=func(x)
dowhile(abs(fx)>1e-6)
f1=func(x1)
if(f1*fx<0)then
x2=x
else
x1=x
endif
x=x2-(x2-x1)/(func(x2)-func(x1))*func(x2)
fx=func(x)
enddo
secant=x
end
realfunctionfunc(x)!
需要求解的函数
realx
func=x**3-2*x**2+7*x+4
end
运行结果:
1.1.2牛顿迭代法
牛顿迭代法基本思路
(1)现任取一个值x1
(2)做一条通过(x1,f(x1))的切线,即以f'(x1)为斜率作直线,直线与x轴的交点为x2,
因为f'(x1)=f(x1)/(x1-x2)
x2=x1-f(x1)/f'(x1)
判断|f(x2)|<ε是否成立,如
果是就找到了这个解,计算完成。
(3)否则,重复步骤
(2),以f'(x1)为斜率做一条通过(x2,f(x2))的切线,直线与x轴的交点为x3,······,直到|f(xn)|<ε,即xn为所得解。
程序编写如下:
realx
integerm
print*,'输入初值'
read*,x
callnewton(x)!
调用牛顿迭代法求解函数
end
subroutinenewton(x)!
牛顿迭代法求解函数
implicitnone
realx,x1
realfunc,dfunc!
对要调用的子程序作说明
integeri,m
i=1
x1=x-func(x)/dfunc(x)
dowhile(abs(x-x1)>1e-6)
print10,i,x1
x=x1
i=i+1
x1=x-func(x)/dfunc(x)
enddo
print20,'x=',x1!
输出计算结果
10format('i=',i4,6x,'x=',f15.7)
20format(a,f15.7)
End
realfunctionfunc(x)!
迭代函数
realx
func=x**3-2*x**2+7*x+4
end
realfunctiondfunc(x)
realx
dfunc=3*x**2-4*x+7
end
运行结果:
1.2求定积分
作业要求:
1、采用函数子程序定义函数f(X);
2、程序选择以下三种方法求定积分:
矩形法、梯形法、辛普生法
3、对于不同的算法分别编写子程序,选择调用,比较不同方法求解的精度。
本题我们用来讨论矩形法、梯形法、辛普生法求定积分的方法。
1.2.1矩形法
矩形法基本思路:
用小矩形面积代替小曲边梯形,矩形面积的求解公式为底×高。
将[a,b]区间分为n个区间,令h=(b-a)/n。
第1个矩形面积:
底=h,高=f(a),也可以用f(a+h)为高,S1=h·f(a)
第i个矩形面积:
底=h,高=f(a+(i-1)·h),也可以用f(a+i·h)为高,Si=h·f(a+(i-1)·h)
程序编写如下:
reala,b,s
integern
realyrectangle
print*,'输入a,b和n的值'
read*,a,b,n
s=rectangle(a,b,n)
print10,a,b,n
print20,s
10format('a=',f5.2,3x,'b=',f5.2,3x,'n=',i4)
20format('s=',f15.8)
End
realfunctionrectangle(a,b,n)
implicitnone
realx,a,b,h,s
integeri,n
realfunc
x=a
h=(b-a)/n
s=0
doi=1,n
s=s+func(x)*h
x=x+h
enddo
rectangle=s
end
realfunctionfunc(x)
realx
func=1+sin(x)
end
运行结果:
n=10时的输出结果
n=100时的输出结果
n=1000时的输出结果
1.2.2梯形法
梯形法基本思路同上,用小梯形面积代替小曲边梯形
第1个梯形面积:
底=h,高=f(a),也可以用f(a+h)为高,S1=h·f(a)
第i个梯形面积:
底=h,高=f(a+(i-1)·h),也可以用f(a+i·h)为高,Si=h·f(a+(i-1)·h)
程序设计如下
reala,b,s
integern
realtrapezia
print*,'输入a,b和n的值'
read*,a,b,n
s=trapezia(a,b,n)
print10,a,b,n
print20,s
10format('a=',f5.2,3x,'b=',f5.2,3x,'n=',i4)
20format('s=',f15.8)
end
realfunctiontrapezia(a,b,n)
implicitnone
realx,a,b,h,s
integeri,n
realfunc
x=a
h=(b-a)/n
s=0
doi=1,n
s=s+(func(x+(i-1)*h)+func(x+i*h))*h/2.0
enddo
trapezia=s
end
realfunctionfunc(x)
realx
func=1+sin(x)
end
运行结果:
n=10时的输出结果
n=100时的输出结果
n=1000时的输出结果
1.2.3辛普生法:
程序编写如下:
reala,b,s
integern
realsinpson
print*,'输入a,b和n的值'
read*,a,b,n
s=sinpson(a,b,n)
print10,a,b,n
print20,s
10format('a=',f5.2,3x,'b=',f5.2,3x,'n=',i4)
20format('s=',f15.8)
end
realfunctionsinpson(a,b,n)
implicitnone
reala,b,h,f2,f4,x
integeri,n
realfunc
h=(b-a)/(2.0*n)
x=a+h
f2=0
f4=func(x)
doi=1,n-1
x=x+h
f2=f2+func(x)
x=x+h
f4=f4+func(x)
enddo
sinpson=(func(a)+func(b)+4.0*f4+2.0*f2)*h/3.0
end
realfunctionfunc(x)
realx
func=1+sin(x)
end
n=10时的输出结果
n=100时的输出结果
n=1000时的输出结果
通过对数据分析我们可以发现,辛普生法最为准确,其次为梯形法,矩形法精度最差
2求解线性方程组
作业要求
用高斯消元法求解线性方程组AX=B的解
其中A为n*n系数矩阵,x为解向量,B为方程组右端维列向量。
要求程序能够求解任意多个未知数的方程组,并附算例及求解结果。
利用Gauss-Jordan法求联立方程组:
有以下联立方程组:
这组等式可以用矩阵方式表示
他们的关系为,c为要求解的未知数。
我们可以先用子程序将矩阵转化为三角矩