牛顿法——求解程序

muzili2008

目的是想建立一个通用的牛顿求解算法的框架，就是说，对于不同的问题，只需要将初值A0，以及问题FX,导数DFX输入，调用牛顿算法，就可以求出结果。


我写的程序如下：有错误，请高手指点一下应该在怎样写？多谢多谢！！

[Fortran] 纯文本查看 复制代码

查看代码

打印?

01	subroutine module_and_algorithm

02	implicit none

03	double precision,external :: cubic_roots,der_cubic_roots

04	integer,parameter :: length=4

05	double precision Co(length)

06	double precision :: x0 =1d0

07	double precision y

08	co=(/3d0,12d0,0d0,1d0/)

09	call newton(x0,cubic_roots,der_cubic_roots,y)

10	print *, "this is the roots near X0",y

11

12	end subroutine

13	subroutine cubic_roots(Co,length,y)

14	implicit none

15	integer length,i

16	double precision x

17	double precision y

18	double precision temp(length)

19	do i=1,length

20	temp(i)=co(i)*x**(i-1)

21

end do

22	y=sum(temp)

23

return

24	end subroutine

25	subroutine der_cubic_roots(Co,length,y)

26	implicit none

27	integer length,i

28	double precision x

29	double precision y

30	double precision temp(length-1)

31	do i=2,length

32	temp(i)=dble(i-1)*co(i)*x**(length-2)

33

end do

34	y=sum(temp)

35

return

36	end subroutine

37	subroutine newton(a,f,df,b)

38	implicit none

39	double precision :: epsilon=1d-6

40	double precision a,b,fb

41	double precision,external :: f,df

42

43	b=a-f(a)/df(a)

44

fb=f(b)

45

46	print *, "i am here (1)",fb

47	do while(abs(fb)>epsilon)

48

a=b

49	b=a-f(a)/df(a)

50

fb=f(b)

51	print *, "i am here (2)", fb

52

end do

53

return

54	end subroutine

55	program main

56	call module_and_algorithm

57

end

muzili2008

解决了  麻烦各位了 。。

li913

初值A0，函数FX,导数DFX，这三者都作为输入参数。
subroutine Newton(a0,fx,dfx,...,...)这样的好处就是代码可不经修改（或少量修改），用于其他项目。