指针之后的问题,包括module、调用其他语言以及GUI,我后面有时间再做吧,感觉对于我来说没啥大的感觉,因为module他其实和c++的class是异曲同工的,没有多少需要用大力气学习的地方,至于调用其他语言,那就更算了,没见过谁Fortran调用c和VB的,感觉没啥用,要调用,也是python调用c++,或者后面看看Python能不能调用Fortran,至于GUI,算了算了,书上的gui画面难看死了,还不如抽时间去学一学Qt,或者是python的canvas、pyqt5和mayavi,再不济学学matlab作图Fortran的作图,就了解一下,知道他还能做gui就行了。
现在主要还是数值计算方法,这个非常重要,虽然Fortran他是又imsl库来解决的,可是你要是不知道数值计算如何进行的,那后面正反演方程你就不会写,更做不出来东西。
一、求解非线性函数:
主要是三种方法:1、二分法(Bisection)2、割线法(Secant)3、牛顿法(Newton,也是最重要的)
1、Bisection:
基本思想:
1、先猜两个值a,b,要让f(a)*f(b)<0,即根据零点定理,再[a,b]之间必有一点c使得f(c)=0
2、取,如果f(c)=0,就找到了正解,如果不等于零,计算f(c),判断f(a)*f(c)以及f(c)*f(b),看哪个小于零,在根据零点定理,重复1、2,直到f(c)=0或者是c的取值达到精度。
来一个
real function func(input)
implicit none
real ::input
func=input-sin(input)
return
end function
program main
implicit none
REAL,PARAMETER::lim=0.0001
REAL i,j,mid
REAL,external::func
! [-1,1]
i=1.0
j=-1.0
do while((func(i)==0).or.(i-j<lim))
if(func(i)==0)then
write(*,*)i
end if
exit
if(func(j)==0)then
write(*,*)j
end if
exit
if(func(i)*func(j)<0)then
mid=(i+j)/2
if(func(mid)*func(i)<0)then
i=mid
else
j=mid
end if
end if
end do
write(*,*)mid
stop "over"
end program
这个程序的输出很简单,因为f(0)=0了,所以输出是这样的:
书上给了一个方程,我们把函数func改一下进行计算:
real function func(input)
implicit none
real ::input
! func=input-sin(input)
func=(input+3)*(input-3)
return
end function
当然这是我的写法,之后我再来把书上给的标准二分法写一遍(这是个模板,改改就能用的那种):
module numerical
implicit none
real,PARAMETER::zero=0.00001
contains
real function bisect(A,B,func)
implicit none
real A,B
real mid
real FA,FB,Fmid
real,external::func
mid=(A+B)/2
Fmid=func(mid)
do while(abs(Fmid)>zero)
FA=func(A)
FB=func(B)
if(FA*FB<0)then
B=mid
mid=(A+B)/2.0
else
A=mid
mid=(A+B)/2.0
end if
Fmid=func(mid)
end do
bisect=mid
return
end function
real function f1(X)
implicit none
real X
f1=X-SIN(X)
end function
end module
program main
use numerical
implicit none
real A,B
real ans
do while(.true.)
write(*,*)"input two number you guess"
read(*,*)A,B
if((f1(A)*f1(B)<0))exit
write(*,*)"wrong area"
end do
! A=-6.0
! B=5.0
ans=bisect(A,B,f1)
write(*,*)ans
stop "over"
end program
咋说呢,他这个程序有问题的,对于这个函数来说,输入的两个数如果是关于原点对称,那就会输出很快,如果输入的上下限不是关于原点对称的,那运行起来就很慢了。
计算起来就很慢。可见,这个二分法在求的时候就不太好用喽。
割线法:
1、先选出两个猜测值a、b
2、画一条通过(a,f(a)),(b,f(b))的一条直线,令这条直线和x轴的交点为c,检查f(c)是否等于0,如果是,就找到了一个解,
3、f(c)!=0的时候,令b、c为两个猜测值a、b,重复上面的步骤。
real function secant(a,b,f)
implicit none
real::a,b,c
real,external::f
real::fa,fb,fc
fa=f(a)
fb=f(b)
c=a-fa*(b-a)/(fb-fa)
fc=f(c)
do while(abs(fc)>zero)
a=b
b=c
fa=f(a)
fb=f(b)
c=a-fa*(b-a)/(fb-fa)
fc=f(c)
end do
secant=c
RETURN
end function
把这段代码加入到上面的module numerical就可以再主程序中调用了,很方便
牛顿法:
1、先做一个猜想值a
2、以作为斜率,结果(a,f(a))作一条直线,令这条直线和x轴的交点为b,检查f(b)=0?
如果是就找到一个解,如果不是则3
3、f(b)!=0的时候,重新令b作为新的猜想值a,重复上述步骤。
咋说呢,主要就是导数求解问题,再就没啥可以说的 了,这个做完了,我都可以用其他的语言来写一个newton,那最后就用python来一个吧(所有的Fortran代码结束)。
real function newtow(a,f,df)
implicit none
real::a
real,external::f
real,external::df
real::b,fb
b=a-f(a)/df(a)
fb=f(b)
do while(abs(fb)>zero)
a=b
b=a-f(a)/df(a)
fb=f(b)
end do
newtow=b
return
end function
real function dfunc(x)
implicit none
real::x
dfunc=COS(x)
return
end function
针对的函数还是把这两个函数加入到module numerical就ok,调用就行
来试一下结果:
ans=nweton(a,func,dfunc)
我没有改输入格式,所以还是输入了两个数,但是真正运算的就是a
这是完整的代码:
module numerical
implicit none
real,PARAMETER::zero=0.0001
contains
real function bisect(A,B,func)
implicit none
real A,B
real mid
real FA,FB,Fmid
real,external::func
mid=(A+B)/2
Fmid=func(mid)
do while(abs(Fmid)>zero)
FA=func(A)
FB=func(B)
if(FA*FB<0)then
B=mid
mid=(A+B)/2.0
else
A=mid
mid=(A+B)/2.0
end if
Fmid=func(mid)
end do
bisect=mid
return
end function
real function secant(a,b,f)
implicit none
real::a,b,c
real,external::f
real::fa,fb,fc
fa=f(a)
fb=f(b)
c=a-fa*(b-a)/(fb-fa)
fc=f(c)
do while(abs(fc)>zero)
a=b
b=c
fa=f(a)
fb=f(b)
c=a-fa*(b-a)/(fb-fa)
fc=f(c)
end do
secant=c
RETURN
end function
real function newtow(a,f,df)
implicit none
real::a
real,external::f
real,external::df
real::b,fb
b=a-f(a)/df(a)
fb=f(b)
do while(abs(fb)>zero)
a=b
b=a-f(a)/df(a)
fb=f(b)
end do
newtow=b
return
end function
real function dfunc(x)
implicit none
real::x
dfunc=COS(x)
return
end function
real function func(x)
implicit none
real::x
func=SIN(x)
return
end function
end module
program main
use numerical
implicit none
real::a,b
real::ans
write(*,*)"input two numbers"
read(*,*)a,b
ans=bisect(a,b,func)!二分法
ans=secant(a,b,func)!割线法
ans=newtow(a,func,dfunc)!牛顿法
write(*,"('X=',F8.4)")ans
stop "over"
end
拿去用吧,记得叫我好心人哈哈,那么到这里就基本上结束了,明天要开始线性代数的使用和学习了。
最后,附上我的python代码和做的图,有点丑,别介意,嘿嘿。
import numpy as np
import matplotlib.pyplot as plt
import math
def df(x)->float:
ans=math.cos(x)
return ans
def func(x)->float:
ans=math.sin(x)
return ans
def newton(a)->float:
b=a-func(a)/df(a)
fb=func(b)
zero=0.00001
while abs(fb)>zero:
a=b
b=a-func(a)/df(a)
fb=func(b)
return b
arr=np.arange(-1,1,0.0001)
arr2=[]
for i in arr:
arr2.append(newton(i))
fig=plt.figure()
plt.plot(arr,arr2)
plt.show()
很丑,哈哈
其实后面我觉得主要还是对牛顿法进行学习,但牛顿法的问题就是求导,得好好求导,确认导数正确。