FORTRAN+计算物理学学习日记（2）

利用Fortran编写数值微分函数，利用向前向后两点式以及五点式求微分。

以上节例题为模板，编写向前向后两点式如下：

!在例一的条件上，计算两点一次插值的微商，编写向前两点式,向后两点式
program mainimplicit noneinteger(4)::i,nreal(4)::x1(7),y1(7),x(100),y(100),p(100),q(100),z(100),hx1(1)=0.0x1(2)=0.5x1(3)=1.0x1(4)=1.5x1(5)=2.0x1(6)=2.5x1(7)=3.0y1(1)=0.0y1(2)=0.4794y1(3)=0.8451y1(4)=0.9975y1(5)=0.9093y1(6)=0.5985y1(7)=0.1411n=6open(unit=101,file="xable.csv",status="new")open(unit=102,file="yable.csv",status="new")open(unit=103,file="pable.csv",status="new")open(unit=104,file="qable.csv",status="new")open(unit=105,file="hable.csv",status="new")h=3.0/100do i=1,100,1x(i)=h*iwrite(101,*) x(i)call interp6(x1,y1,x(i),n,y(i))write(102,*)y(i)end dodo i=1,99,1p(i)=(y(i+1)-y(i))/hwrite(103,*) p(i)end doq(1)=(y(1)-0.00)/hwrite(104,*)q(1)   do i=2,100,1q(i)=(y(i)-y(i-1))/hwrite(104,*) q(i)end do   close(101)close(102)close(103)close(104)end program mainsubroutine interpolyb(x1,x,n,yy)
integer k,j
integer,intent(in)::n
real,intent(in)::x1(n+1)
real,intent(in)::x
real,intent(out)::yy(n+1)
yy(1)=1;
do k=2,n+1yy(1)=yy(1)*(x-x1(k))/(x1(1)-x1(k))
end do
yy(n+1)=1;
do k=1,nyy(n+1)=yy(n+1)*(x-x1(k))/(x1(n+1)-x1(k));
end do
do j=2,nyy(j)=1do k=1,j-1yy(j)=yy(j)*(x-x1(k))/(x1(j)-x1(k))end dodo k=j+1,n+1 yy(j)=yy(j)*(x-x1(k))/(x1(j)-x1(k)); end do
end do
return end subroutine interpolybsubroutine interp6(x1,y1,x,n,y)
real,intent(out)::y
real yy(n+1)
integer::i
integer,intent(in)::n
real,intent(in)::x1(n+1)
real,intent(in)::y1(n+1)
real,intent(in)::x
y=0
call interpolyb(x1,x,n,yy)
do i=1,n+1y=y+yy(i)*y1(i)
end do
end subroutine interp6

绘制如下图片：

五点式代码及绘制图像：

!在例一的条件上，改动部分函数，编写五点式方程
program mainimplicit noneinteger(4)::i,nreal(4)::x1(7),y1(7),x(103),y(103),z(103),hx1(1)=0.0x1(2)=0.5x1(3)=1.0x1(4)=1.5x1(5)=2.0x1(6)=2.5x1(7)=3.0y1(1)=0.0y1(2)=0.4794y1(3)=0.8451y1(4)=0.9975y1(5)=0.9093y1(6)=0.5985y1(7)=0.1411n=6open(unit=101,file="xxable.csv",status="new")open(unit=102,file="yyable.csv",status="new")open(unit=103,file="hable.csv",status="new")h=3.0/100do i=1,103,1x(i)=h*i-2*hwrite(101,*) x(i)call interp6(x1,y1,x(i),n,y(i))write(102,*)y(i)end doclose(101)close(102)!编写五点式do i=3,101,1z(i)=(y(i-2)-8*y(i-1)+8*y(i+1)-y(i+2))/(12*h)write(103,*)z(i)end doclose(103)end program mainsubroutine interpolyb(x1,x,n,yy)
integer k,j
integer,intent(in)::n
real,intent(in)::x1(n+1)
real,intent(in)::x
real,intent(out)::yy(n+1)
yy(1)=1;
do k=2,n+1yy(1)=yy(1)*(x-x1(k))/(x1(1)-x1(k))
end do
yy(n+1)=1;
do k=1,nyy(n+1)=yy(n+1)*(x-x1(k))/(x1(n+1)-x1(k));
end do
do j=2,nyy(j)=1do k=1,j-1yy(j)=yy(j)*(x-x1(k))/(x1(j)-x1(k))end dodo k=j+1,n+1 yy(j)=yy(j)*(x-x1(k))/(x1(j)-x1(k)); end do
end do
return end subroutine interpolybsubroutine interp6(x1,y1,x,n,y)
real,intent(out)::y
real yy(n+1)
integer::i
integer,intent(in)::n
real,intent(in)::x1(n+1)
real,intent(in)::y1(n+1)
real,intent(in)::x
y=0
call interpolyb(x1,x,n,yy)
do i=1,n+1y=y+yy(i)*y1(i)
end do
end subroutine interp6

FORTRAN+计算物理学学习日记（2）相关推荐

FORTRAN+计算物理学学习日记（1）
第一周:结合李录的计算物理学学习FORTRAN语言,这周的任务是插值函数,大致编写了四个小时,编写了一个双层循环的插值函数,进行了六次插值计算例题. 例题如下: 编写代码如下: program mai ...
FORTRAN+计算物理学学习日记（5）
2.1常微分方程的简单数值解法本节编写了四种简单的数值方法去求解常微分方程的初始问题,包括 Euler 方法.Taylor 级数法.后向 Euler 方法和梯形公式. 注意点:定义格式或者数组定义 ...
FORTRAN+计算物理学学习日记（6）
2.2Runge-Kutta 方法求解常微分方程 "直接利用 Taylor 级数展开提高算法的阶数有许多困难,特别是要确定函数 f (x, y) 的导数,这在数值计算中是非常不方便的．为了 ...
FORTRAN+计算物理学学习日记（7）
2.3多步法求解常微分方程 program mainimplicit noneinteger(8)::n,m,i,kreal(8),external::fa,fbreal(8)::h,x(200),y ...
FORTRAN+计算物理学学习日记（4）
1.5基本数学运算中的求根方法一:区间对分法求根书中例题及编写代码如下 !!利用区间对分法求根 program mainimplicit nonereal(8)::a,b,x,t,ya=2b=3x ...
FORTRAN+计算物理学学习日记（8）
第三章边值问题和本征值问题 3.1numerov算法例题 program mainimplicit noneinteger(8)::i,nreal(8)::x(600),y(600),a,h,pi ...
Java学习日记-Day01
Java学习日记-Day01 Java语言概述比特(byte)与字节内存 Java基础知识图解人机交互方式常用的DOS命令常用快捷键计算机编程语言介绍第一代语言第二代语言第三代语言 ...
深度学习日记 2 - 概率论与信息论基础
深度学习日记 2 - 概率论与信息论基础: 1.随机变量(random variable):是可以随机地取不同值的变量.我们通常用打印机体的小写字母来表示随机变量本身,而用脚本字体中的小写字母来表示 ...
GPU（CUDA）学习日记（十一）------ 深入理解CUDA线程层次以及关于设置线程数的思考
GPU(CUDA)学习日记(十一)------ 深入理解CUDA线程层次以及关于设置线程数的思考标签: cuda存储线程结构网格 2012-12-07 16:30 6298人阅读评论(4)收藏举 ...

FORTRAN+计算物理学学习日记（2）

FORTRAN+计算物理学学习日记（2）相关推荐

最新文章

热门文章