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

常微分方程的边值问题和本征值问题一、 问题描述利用搜索法和弦割法，得到该常微分方程的本征值，再利用打靶法计算多个本征值。二、 解决方法（一）搜索法1.先随便猜测k的一个试验值,程序中令k=12.由Numerov算法根据本题的条件，kn+1=kn=kn-1=k,s=0,得到yn+2,yn+1,yn间的迭代公式令con=(k*h)2/12yn+2=2*(1-5*con)*yn+1/(1+con)-yn3自己给定的初始条件，然后利用公式得到边界值(1)4.然后以小的步长dk增加k值，这里令dk=1,每当(1)改变符号时，就将步长减半后倒退回来重复5.当步长小于所要求的容许误差时终止程序，此时的k值即为所求。（二） 弦割法1. 随便猜测两个k值，这里令k0=1,k1=22. 自己给定的初始条件，对两个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 returnend subroutine addprogram 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*b10)then !若果y1变号 k=k-dk !后退 dk=dk/2.0 !步长减半 end if end do write(*,*)k !写出求得的本征值end2. 弦割法 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点的位置做准备 returnend 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)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 dowrite(*,*)k2 !得到本征值k end program zy33. 打靶法得多个本征值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 returnend subroutine addprogram 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*b10)then !若果y1变号 k=k-dk !后退 dk=dk/2.0 !步长减半 end if end do write(*,*)k !写出求得的本征值 end doend4. 选取一个本征值，看函数图像（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.001open(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 doclose(10) end四、 程序结果1. 搜索法2. 弦割法3. 打靶法4. k=1.5708时的函数五、 实验中出现的问题在利用打靶法找本征值的过程中，如果步长太小，会出现如下情况而当笔者将3代入k值用搜索法的时候，得到笔者认为出线该情况的可能是循环程序彼此间出现了问题，而当调整了步长之后，冲突减小，得到比较正常的特征值