计算物理实验报告常微分方程.doc

常微分方程的边值问题和本征值问题

一、 问题描述

利用搜索法和弦割法，得到该常微分方程的本征值，再利用打靶法计算多个本征值。

二、 解决方法
（一）搜索法
1.先随便猜测k的一个试验值,程序中令k=1
2.由Numerov算法

根据本题的条件，kn+1=kn=kn-1=k,s=0,得到yn+2,yn+1,yn间的迭代公式
令con=(k*h)^2/12
yn+2=2*(1-5*con)*yn+1/(1+con)-yn
3自己给定φ的初始条件，然后利用公式得到边界值φ(1)
4.然后以小的步长dk增加k值，这里令dk=1,每当φ(1)改变符号时，就将步长减半后倒退回来重复
5.当步长小于所要求的容许误差时终止程序，此时的k值即为所求。

（二） 弦割法
1. 随便猜测两个k值，这里令k0=1,k1=2
2. 自己给定φ的初始条件，对两个k值分别利用上述公式进行迭代，得到边界值y1(1)和y2(1)。
3. 比较y1(1)和y2(1)的绝对值大小。若绝对值大，说明对应的k值距离本征值距离较远。
4. 将（k0+k1）/2赋给k2,边界值绝对值小的对应的k值保持不变，边界值绝对值大的对应k值重新定位k2的值。
5. 重复进行实验，当y1(1)和y(2)的差的绝对值小于容许误差时终止程序。此时k1的值即为所求。

当搜索法和弦割法大致求出了一个本征值后，利用打靶法，调整k值再度进行搜索，得到多个本征值，绘出其中一个本征值对应的函数图像，观察其性质。

三、 程序实现
1. 搜索法
subroutine add(t,y0,y1) !利用子程序表示函数值的迭代
implicit none
real(8)::t,h,con,y0,y1,y2
integer::i,n
n=10000
h=1.0/n
con=(t*h)**2/12
do i=1,n-1
y2=2*(1-5*con)*y1/(1+con)-y0 !利用Numerov算法，得到迭代公式
y0=y1 !向前迭代
y1=y2
end do
return
end subroutine add
program zy3
implicit none
real(8)::diffk,dk,yold,k,b0,b1
integer::s
b0=0.01 !取初始值，根据题目条件，令y0=y1,来保证x=0的位置导数为0
b1=0.01
s=1
k=s !给定一个猜测的k值,此为搜索的初值
dk=1 !给定步长
diffk=0.0000001 !给定步长最后达到的误差范围
call add(k,b0,b1)
yold=b1 !通过运行子程序，得到由初始值积到x=1时的不为0的函数值
do while(abs(dk)>diffk) !开始搜索
k=k+dk !在k中走一步
b0=0.01
b1=0.01
call add(k,b0,b1)
if(yold*b1<0)then !若果y1变号
k=k-dk !后退
dk=dk/2.0 !步长减半
end if
end do
write(*,*)k !写出求得的本征值
end
2. 弦割法
subroutine add(t,b) !利用子程序表示函数值的迭代
implicit none
real(8)::t,h,con,y0,y1,y2,b
integer::i,n
b=0
y0=0.01
y1=0.01
n=10000
h=1.0/n
con=(t*h)**2/12
do i=1,n-1
y2=2*(1-5*con)*y1/(1+con)-y0 !利用Numerov算法，得到迭代公式
y0=y1 !向前迭代
y1=y2
end do
b=abs(y1) !得到x=1处函数值的绝对值，为确定k2点的位置做准备
return
end subroutine add
program zy3
real(8)::a,k0,k1,k2,dk,m1,m2,dm
integer::i,n
k0=1 !给两个启动值
k1=2
k2=0
dm=0.00000001 !表示k0和k1对应的函数值相等时允许的误差
m1=0 !此值表示k值取k0时，x=1处函数值的绝对值
m2=0 !此值表示k值取k1时，x=1处函数值的绝对值
do while (.true.)
call add(k0,m1)
call add(k1,m2) !运行子程序，分别得到两个k值对应的x=1处的函数值的绝对值
if(abs(m1-m2)<dm)exit !当k0k1对应的绝对值近似相等时退出循环
if(m1>m2)then !如果k0对应的函数值绝对值较大
k2=(k1+k0)/2.0 !k2点取在k1和k0的平均值
k0=k2 !当k0对应的函数值绝对值较大时，表示其离本征值较远，而将其舍弃不用，赋k2值
k1=k1 !此时k1距离本征值较近，不用变
else
k2=(k1+k0)/2.0 !反之，舍去k1，k0不变
k0=k0
k1=k2
end if

end do
write(*,*)k2 !得到本征值k
end program zy3
3. 打靶法得多个本征值
subroutine add(t,y0,y1) !利用子程序表示函数值的迭代
implicit none
real(8)::t,h,con,y0,y1,y2
integer::i,n
n=10000
h=1.0/n
con=(t*h)**2/12
do i=1,n-1
y2=2*(1-5*con)*y1/(1+con)-y0 !利用Numerov算法，得到迭代公式
y0=y1 !向前积分
y1=y2
end do
return
end subroutine add
program zy3
implicit none
real(8)::diffk,dk,yold,k,b0,b1
integer::s
b0=0.01 !取初始值，根据题目条件，令y0=y1,来保证x=0的位置导数为0
b1=0.01
do s=1,100,2 !改变k的初值
k=s
dk=1 !给定步长
diffk=0.0000001 !给定步长最后达到的误差范围
call add(k,b0,b1)
yold=b1 !通过运行子程序，得到由初始值积到x=1时的不为0的函数值
do while(abs(dk)>diffk) !开始搜索
k=k+dk !在k中走一步
b0=0.01
b1=0.01
call add(k,b0,b1)
if(yold*b1<0)then !若果y1变号
k=k-dk !后退
dk=dk/2.0 !步长减半
end if
end do
write(*,*)k !写出求得的本征值
end do
end
4. 选取一个本征值，看函数图像（k=1.5708）
program zy3
implicit none
real(8)::diffk,dk,yold,k,y0,y1,y2,h,con
integer::i,n
k=1.5708
y0=0.001
y1=0.001
open(unit=10,file='p.txt')
n=1000
h=10.0/n
con=(k*h)**2/12
do i=1,n-1
write(10,*)y1,i*h
y2=2*(1-5*con)*y1/(1+con)-y0
y0=y1
y1=y2
end do
close(10)
end
四、 程序结果
1. 搜索法

2. 弦割法

3. 打靶法

4. k=1.5708时的函数

五、 实验中出现的问题
在利用打靶法找本征值的过程中，如果步长太小，会出现如下情况

而当笔者将3代入k值用搜索法的时候，得到

笔者认为出线该情况的可能是循环程序彼此间出现了问题，而当调整了步长之后，冲突减小，得到比较正常的特征值