计算物理实验报告 常微分方程 下载本文

内容发布更新时间 : 2024/5/18 23:32:18星期一 下面是文章的全部内容请认真阅读。

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

一、问题描述

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

二、解决方法 (一)搜索法

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 b1=0.01 s=1

k=s dk=1 diffk=0.0000001 call add(k,b0,b1)

yold=b1 do while(abs(dk)>diffk) k=k+dk b0=0.01 b1=0.01

call add(k,b0,b1)

if(yold*b1<0)then 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

取初始值,根据题目条件,令y0=y1,来保证x=0的位置

导数为0

!给定一个猜测的k值,此为搜索的初值 !给定步长

!给定步长最后达到的误差范围 !通过运行子程序,得到由初始值积到x=1时的不为0

的函数值

!开始搜索 !在k中走一步 !若果y1变号 !后退 !步长减半 !写出求得的本征值

!利用子程序表示函数值的迭代 ! 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)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