常微分方程边值问题的试射法毕业论文.doc

南 京 师 范 大 学

毕 业 设 计（论 文）
（ 2011 届）

题 目： 常微分方程边值问题的试射法
学 院： 数学科学学院
专 业： 信息与计算科学
姓 名：
学 号：
指导教师：

目录
目录 1
常微分方程边值问题的试射法 2
摘要 2
关键字 2
Abstract 3
Keyword 3
1.引言 4
2．二阶常微分方程第一边值问题 4
2.1二阶常微分方程边值问题的相关定理 4
2.2二阶常微分方程第一边值问题的算法 6
2.3使用上述算法求解 10
3.二阶常微分方程第二，第三边值问题 12
3.1二阶常微分方程第二边值问题算法思想和步骤 12
3.2二阶常微分方程第三边值问题算法思想和步骤 14
4.其余常用算法的概述 15
4.1 Netwon法的概述 15
4.2二分法的概述 15
4.3三种方法的比较 16
5算法的改进 16
6.致谢词 17
7.参考文献 17
附录1 二阶常微分方程第一边值问题的试射法（割线法）程序（c++） 19
附录2 英文文献翻译 22

常微分方程边值问题的试射法
张晓飞
06级3班 信息与计算科学专业 06070312

摘要
本论文给出了二阶常微分方程边值问题的相关概念及描述，重点讨论了二阶常微分方程第一边值问题。并给出了二阶常微分方程第一边值问题解的存在唯一性定理，以及相关的定理证明。利用化边值问题为初值问题的基本思想，给出了二阶常微分方程第一，第二，第三边值问题试射法及打靶法的基本思想和相关的算法步骤。给出了二阶常微分方程第一边值问题试射法的实例和程序，最后并对试射法进行了些许的改进。
关键字
常微分方程边值问题 试射法 割线法 牛顿法

Abstract
This paper gives the description about the related concepts of second order ordinary differential equations, focusing on the first boundary value problem of second order ordinary differential equations. Then gives the first order ordinary differential equations boundary value problem and uniqueness theorem and proved the theorems. And I give the basic idea and the related steps of the algorithm of shooting method , which work out the first, second and third boundary value about the second order ordinary differential equations, using the basic idea about using initial value problem instead of boundary value problem. I gives a instances and a program about the first boundary value problem of second order ordinary differential equations. Finally, I do a little improvement about the shooting method.
Keyword
Ordinary Differential Equations shooting method Secant method Newton

1.引言
常微分方程边值问题在应用科学与工程技术中是经常遇到的，由于边值问题解的存在唯一性比初值问题复杂多，至今没有很有效的解决方法，因此在解决边值问题无论在理论还是在实际计算中都比初值问题麻烦。采用较好的计算格式来进行计算，不仅可以提高解决问题的效率，更能提高应用计算结果的精度，从而使之更能广泛地应用于科学工程领域。而常微分方程的试射法就是被广泛采用的数学方法。
2．二阶常微分方程第一边值问题
二阶常微分方程第一边值问题
（1）
将问题（1）转化为初值问题
（2）
的解在的值满足或者其中为允许的误差界。这样，我们就把作为边值问题（1）的近似解。显然试射法最终将边值问题化为初值问题，再利用相应的方法求解初值问题。在此过程中最为重要的是找出最为合适的，尽可能多的减少计算步骤。
2.1二阶常微分方程边值问题的相关定理
2.1.1二阶线性常微分方程解的存在唯一性定理
（3）
满足
（I），，在上连续；
（II）>0,对。
则边值问题(3)的解存在唯一。
由于非齐次问题的解等于其相应的齐次问题的通解加上其本身的一个特解之和。对于二阶线性边值问题（3）的解可以用以下非齐次线性初值问题（4）的解和齐次线性初值问题的解来表述出来。
（4）
（5）
证明：设是非齐次线性初值问题（4）的解，是齐次线性问题（5）的解，设

就是二阶非线性边值问题的唯一解。
2.1.2二阶常微分方程解的存在唯一性定理
设方程（1）的函数及，在区间内连续，且
（I）
（II）在内有界，即存在常数M，使得

则边值问题（1）的解存在唯一。
2.2二阶常微分方程第一边值问题的算法
2.2.1二届常微分方程第一边值问题的算法思想
只要问题（2）的解在的值满足或者，为允许误差界。这样，我们把作为边值问题（1）的近似解。为此，可以采取逐次逼近法来实现。
假设为边值问题的解，我们估计的值为后，解决初值问题

这样得到的解为，并计算得到。一般，若或者，则把作为边值问题（1）的近似解；否则必须调整，例如取（），则在（2）中令，再求解此初值问题。设计算得它的解为，，若或，则作为边值问题（1）的近似解。否则再修改。如此重复计算，直至或者时，便以作为边值问题（1）的近似解
y
x
O
a
b

图1
参见图1上面的积分曲线的过大，下面的积分曲线的过小。
对于解高阶的微分方程初值问题的数值法可以使用经典四阶Runge-Kutta方法求解。则问题关键是如何确定参数。当，解得初值问题的解为。自然，我们希望

从而，确定的问题可以归结为求方程
（6）
的近似跟。（6）是一个非线性方程，可以使用二分法，割线法，Netwon法来求解，本文只采用割线法和二分法。
2. 2.2二阶常微分方程的第一边值问题算法过程
2.2.2.1经典四阶Runge Kutta法求解高阶微分方程初值问题算法
求解高阶微分方程初值问题，即将高阶化为一阶常微分方程组的初值问题

输入：端点a，b；方程个数m；整数N；初值，，·····，。
输出：解在t的N+1个等距点的近似。
Step1. ；
。
Step2. 对j=1，2，·····，m
。
Step3. 输出。
Step4. 对i=1，2，·····，N做Step5~11.
Step5. 对j=1，2，···，m
。
Step6. 对j=1，2，···，m
。
Step7. 对j=1，2，···，m

Step8. 对j=1，2，···，m

Step9. 对j=1，2，···，m

Step10. 。
Step11. 输出。
Step12. 停机。
2.2.2.2二阶常微分边值问题的试射法（割线法）
选取割线法来解问题（6）我们需要选取初始近似的和，由公式

生成序列｛｝。按照割线法求得，直到为止，其中为允许的误差界。
因此，解边值问题（1）的打靶法（由割线法确定参数，由经典四阶Runge-Kutta方法求解）的计算步骤如下：
令。给定初始值，误差容限，最大迭代次数。
Step1. 取，解初值问题（2）得解在的近似值：
，·····，，（=）
Step2. 若则输出
，········，，
作为初值问题（1）的解在，·····，，的近似值；停机。否则转Step3。
Step3. 令，解初值问题（2）得解在的近似值：
，·······，，（= ）。
Step4. 若，则输出
，·······，，
作为初值问题（1）的解在，·····，，的近似值；停机。否则转Step5。
Step5. 对于k=2，3，·····，m做 Step6~8
Step6. 由割线法公式计算；
Step7. 令，解初值问题（2）得解在（=1，·····N）的（近似）值：
，·······，，（=）
Step8. 若，则输出
，·····，，
作为初值问题（1）的解在，·····，，的近似值；停机。
Step9. 输出（‘Method failed’）；
停机。
2.3使用上述算法求解
2.3.1当所取初始值为正数时
初始输入如图1：

图1

结果如表1：
i
x(i)
数值解
准确解
误差
0
0
-0.3
-0.3
0
1
pi/20
-0.311949732
-0.311949949
2.167E-07
2
pi/10
-0.316218321
-0.316218654
3.338E-07
3
3pi/20
-0.31270064
-0.312701007
3.675E-07
4
pi/5
-0.301483289
-0.301483624
3.35E-07
5
pi/4
-0.282842458
-0.282842713
2.55E-07
6
3pi/10
-0.257237127
-0.257237275
1.487E-07
7
7pi/20
-0.225297763
-0.225297802
3.98E-08
8
2pi/5
-0.187810794
-0.18781075
-4.41E-08
9
9pi/20
-0.145699244
-0.145699174
-7.05E-08
10
pi/2
-0.1
-0.1
0

表1
此时试射法中t的取值为：-0.09999343444
一共计算3次。
精确解与近似解的图像如图2：
其中黑色曲线表示精确解，蓝色曲线表示近似解。
图2

2.3.2当所取初始值为负数时
初始输入如图3：

图3
结果如表2：

i
x(i)
数值解
准确解
误差
0
0
-0.3
-0.3
0
1
pi/20
-0.311949732
-0.311949949
2.167E-07
2
pi/10
-0.316218321
-0.316218654
3.338E-07
3
3pi/20
-0.31270064
-0.312701007
3.675E-07
4
pi/5
-0.301483289
-0.301483624
3.35E-07
5
pi/4
-0.282842458
-0.282842713
2.55E-07
6
3pi/10
-0.257237127
-0.257237275
1.487E-07
7
7pi/20
-0.225297763
-0.225297802
3.98E-08
8
2pi/5
-0.187810794
-0.18781075
-4.41E-08
9
9pi/20
-0.145699244
-0.145699174
-7.05E-08
10
pi/2
-0.1
-0.1
0

表2
此时试射法中t的取值为：-0.09999343444
一共计算3次。
精确解与近似解的图像如图4：
其中黑色曲线表示精确解，蓝色曲线表示近似解。

图4
3.二阶常微分方程第二，第三边值问题
3.1二阶常微分方程第二边值问题算法思想和步骤
（7）
将问题（7）转换为问题出边值问题（8）
（8）
的解在的值满足或者其中为允许的误差界。这样，我们就把作为边值问题（7）的近似解。具体求解步骤如下：
令。给定初始值，误差容限，最大迭代次数。
1.（1）取解初值问题（8）得到在，处的近似解
（2）若则输出：
，········，，
作为初值问题（1）的解在，·····，，的近似值；停机。
2. （1）令，解初值问题（2）得解在的近似值
（2）若则输出：
，·······，，
作为初值问题（1）的解在，·····，，的近似值；停机。
3.对做
（1）用割线法计算得
（2）令，解初值问题（2）得解在（=1，·····N）的（近似）值
（3）若则输出：
，·····，，
作为初值问题（1）的解在，·····，，的近似值；停机。
3.2二阶常微分方程第三边值问题算法思想和步骤
（9）
将问题（9）转换为问题出边值问题（10）设
（10）
的解在的值满足或者其中为允许的误差界。这样，我们就把作为边值问题（10）的近似解。
其试射法基本过程如下：
（1） 选取参数，令。设，则由可以确定，从而得到初值问题
（11）
求初值问题（11）的解
（2）把代入的左端得到
。
若，则为所求边值问题的近似解；否则，则令

设（用代替）解初值问题（11），求得解后再代入的左端得
。
（3）从，，，出发，由割线法

产生，以代替解初值问题（10）确定，如此继续执行下去，得到序列,直到

为止，为误差容限。
4.其余常用算法的概述
4.1 Netwon法的概述
求解方程的一个根，利用牛顿法，我们只需选取一个，由公式

生成序列，直到为止，为允许误差界。
具体算法步骤中，只需在求出的同时，求出。
4.2二分法的概述
由函数零点定理可知，只要连续函数，满足，则在中一定一个点使得函数值为0。任取两点，使得，再由公式

生成序列，直到为止，为允许误差界。
在选取二分区间时可由下述定理确定具体区间
定理：对于初值问题中及，在区间内连续，并且
（I）
（II）存在常数M，使得

则初值问题的解关于严格单调增加。
4.3三种方法的比较
割线法：优点：与牛顿法相比，割线法的每一步只要计算一次函数值，而牛顿法还需要计算导数值。因此在导数计算比较费事或不可能的情况下，割线法则显出其优点。
缺点：割线法的收敛阶为1.618阶，一般情况下比牛顿法的收敛速度慢。
牛顿法：优点：对于一般问题牛顿法的收敛阶为2阶，其收敛速度快。
缺点：每一步除了要计算一次函数值外，还要计算导数值。
二分法：优点：直观，简单，安全可靠，迭代速度也快。
缺点：近似解的误差下降速度不快，而且对于未知的方程，找到适当的和使得，不是那么容易。
5算法的改进
从上述内容可以看到试射法的关键在于如何高效确定出数列，使得能在尽量短的时间中，求出方程的近似解。在确定时，关键分为两步：
（1） 找出一个适当的；
（2） 找出一个高效的公式，求出。
由于二分法的近似解误差下降速度不快，所以可以先用二分法找到一个，将此作为其他方法的。
或者可将该问题转化为无约束最优化问题

利用解无约束最优化的问题的方法，来求解。
6.致谢词
在这次毕业设计的整个过程中，虽然遇到了一些技术上的难题但在与老师、同学和网友的交流下都及时的解决了问题。其中魏虹老师认真负责的工作态度，严谨的治学精神和深厚的理论水平都使我受益匪浅。她无论在理论上还是在实践中，都给与我很大的帮助，使我得到不少的提高，这对于我以后的工作和学习都有一种巨大的帮助，感谢她耐心的辅导。
我经过几个月的收集资料、阅读文献、写作论文，终于顺利完成了论文。我要感谢很多人，首先感谢学校和老师们，不仅为我们提供了良好的设施和学习环境，而且教给我们很多知识和做人的道理。还要感谢论文指导老师魏虹老师，从课题的选择到最终完成，魏老师都始终给予我细心的指导和不懈的支持。魏老师认真负责的工作态度，严谨的治学精神和深厚的理论水平都使我受益匪浅。她无论在理论上还是在实践中，都给与我很大的帮助，使我得到不少的提高，这对于我以后的工作和学习都有一种巨大的帮助，感谢她耐心的辅导，在此谨向魏老师致以诚挚的谢意和崇高的敬意。还要感谢辅导员王海康和喻华杰老师，多年来一直照顾我们的学习和生活。还要感谢在一起愉快地度过大学4年学习生活的各位同学，特别是我的室友们，正是由于你们的帮助和支持，我才能克服一个又一个的困难和疑惑，直至本文的顺利完成。同时感谢我的同组成员，给予了我不少的帮助。最后我要感谢含辛茹苦培养我长大的父母，有你们才有今天的我。谢谢你们！
通过此次论文，我学到很多知识。通过收集资料，增强了我的自学和动手能力，通过论文写作，不仅提高了我撰写科技论文和技术报告的能力，而且培养了我的创新意识和严谨 的科学作风。
即将挥别我的学校、老师、同学，心里只有千万不舍，但是对未来的路，我充满了信心。最后，感谢在大学期间认识我和我认识的所有人，有你们伴随，才使得我大学的生活如此丰富多彩，绚丽多姿！
7.参考文献
【1】余德浩，汤华中.微分方程数值解法.科学出版社.2005
【2】胡健伟，汤怀民.微分方程数值解法.科学出版社.1999
【3】吴微.非线性抛物型方程广义差分法的误差估计.计算数学.1987
【4】李荣华，冯果忱.微分方程数值解法.高等教育出版社.1988
【5】李荣华，冯果忱.微分方程数值解法(第三版).高等教育出版社.1996
【6】李庆阳，王能超，易大义.数值分析. 北华中科技大学出版社.1986
【7】姜礼尚，陈亚浙，刘西垣，易法槐.数学物理方程讲义.高等教育出版社.2007
【8】周顺兴.解抛物型偏微分方程的高精度差分格式.计算数学.1982
【9】Stig Larsson, Vidar Thomée..Partial differential equations with numerical methods .
科学出版社.2006
【10】谭浩强等.C程序设计.清华大学出版社.1999
【11】钱能.C++程序设计教程（第二版）.清华大学出版社.2005
【12】Richard L.Burden，J.Douglas Faires，Numerical Analysis(Seventh Edi-
tion)(影印版). 高等教育出版社.2002
【13】矢岛信男，野木达夫.发展方程的数值解析.东京岩波书店.1977

附录1 二阶常微分方程第一边值问题的试射法（割线法）程序（c++）
#include<iostream>
#include<cmath>
#include<iomanip>
#include<fstream>
using namespace std;
double f(double a,double b,double c)
{
return cos(a)+2*b+c;
}
void fourRungeKutta0(double a,double b,double N,double c[2])
{
double h,t;
double w[2];
double k[4][2];
h=(b-a)/N;
t=a;
w[0]=c[0];
w[1]=c[1];
ofstream fout("a.txt");
fout<<setprecision(10)<<w[0]<<endl;
for (int i=0;i<N;i++)
{
k[0][0]=h*w[1];
k[0][1]=h*f(t,w[0],w[1]);
k[1][0]=h*(w[1]+0.5*k[0][1]);
k[1][1]=h*f(t+0.5*h,w[0]+0.5*k[0][0],w[1]+0.5*k[0][1]);
k[2][0]=h*(w[1]+0.5*k[1][1]);
k[2][1]=h*f(t+0.5*h,w[0]+0.5*k[1][0],w[1]+0.5*k[1][1]);
k[3][0]=h*(w[1]+k[2][1]);
k[3][1]=h*f(t+h,w[0]+k[2][0],w[1]+k[2][1]);
w[0]=w[0]+(k[0][0]+2*k[1][0]+2*k[2][0]+k[3][0])/6.0;
w[1]=w[1]+(k[0][1]+2*k[1][1]+2*k[2][1]+k[3][1])/6.0;
t=a+(i+1)*h;
fout<<w[0]<<endl;
// cout<<"下一步"<<endl;
}

}
double fourRungeKutta(double a,double b,double N,double c[2])
{
double h,t;
double w[2];
double k[4][2];
h=(b-a)/N;
t=a;
w[0]=c[0];
w[1]=c[1];
// cout<<"("<<t<<","<<w[0]<<","<<w[1]<<")"<<endl;
for (int i=0;i<N;i++)
{
k[0][0]=h*w[1];
k[0][1]=h*f(t,w[0],w[1]);
k[1][0]=h*(w[1]+0.5*k[0][1]);
k[1][1]=h*f(t+0.5*h,w[0]+0.5*k[0][0],w[1]+0.5*k[0][1]);
k[2][0]=h*(w[1]+0.5*k[1][1]);
k[2][1]=h*f(t+0.5*h,w[0]+0.5*k[1][0],w[1]+0.5*k[1][1]);
k[3][0]=h*(w[1]+k[2][1]);
k[3][1]=h*f(t+h,w[0]+k[2][0],w[1]+k[2][1]);
w[0]=w[0]+(k[0][0]+2*k[1][0]+2*k[2][0]+k[3][0])/6.0;
w[1]=w[1]+(k[0][1]+2*k[1][1]+2*k[2][1]+k[3][1])/6.0;
t=a+(i+1)*h;
// cout<<"("<<t<<","<<w[0]<<","<<w[1]<<")"<<endl;
// cout<<"下一步"<<endl;
}
return w[0];
}
void main()
{
double a,b,tol,m,z,y;
double N;
double c[2],t[100],d[100];
// double w[2];//d[100]中放置每次取一个t后得到的最后解w[0][N-1]
cout<<"输入a的值：a=";
cin>>a;
cout<<"输入b的值：b=";
cin>>b;
cout<<"输入初始t[0]：t[0]=";
cin>>t[0];
cout<<"输入最大迭代次数m：m=";
cin>>m;
cout<<"输入误差容限tol：tol=";
cin>>tol;
cout<<"输入分割数N：N=";
cin>>N;
cout<<"输入精确解在a点的值y：y=";
cin>>y;
cout<<"输入精确解在b点的值z：z=";
cin>>z;
c[0]=y;c[1]=t[0];
cout<<setprecision(10)<<a<<endl;
cout<<"t[0]="<<t[0]<<endl;
d[0]=fourRungeKutta(a,b,N,c);
cout<<"d[0]="<<d[0]<<endl;
if(fabs(d[0]-z)<tol)
{
fourRungeKutta0(a,b,N,c);
}
else
{
t[1]=z/d[0]*t[0];
c[1]=t[1];
cout<<"t[1]="<<t[1]<<endl;
d[1]=fourRungeKutta(a,b,N,c);
cout<<"d[1]="<<d[1]<<endl;
}
if(fabs(d[1]-z)<tol)
{
fourRungeKutta0(a,b,N,c);
}
else
{
for(int k=2;k<=m;k++)
{
t[k]=t[k-1]-(d[1]-z)*(t[k-1]-t[k-2])/(d[1]-d[0]);
c[1]=t[k];
d[k]=fourRungeKutta(a,b,N,c);
cout<<"t[k]="<<t[k]<<endl;
cout<<"d[k]="<<d[k]<<endl;
if(fabs(d[k]-z)<tol)
{
cout<<"输出了："<<endl;
fourRungeKutta0(a,b,N,c);
break;
}

}
cout<<k;
}
}
附录2 英文文献翻译
英文：
8.4 A Maximum Principle
We now consider the generalization of the mixed initial-boundary value problem of Sect. 8.2 which allows a source term and inhomogeneous boundary conditions, i.e., to find such that
(8.34)
Where is a bounded domain in and is a finite interval in time. In order to show a maximum principle for this problem it is convenient to introduce the parabolic boundary of as the set ，i.e., the boundary of minus the interior of the top part of this boundary, .
Theorem 8.6. Let u be smooth and assume that in .Then u attains its maximum on the parabolic boundary.
Proof. If this were not true, then the maximum would be attained either at an interior point of or at a point of, i.e., at a point , and we would have
.
In such a case, for sufficiently small, the function

would also take its maximum at a point in, since, for small,

By our assumption we have since that
(8.35)
On the other hand, at the point, where w takes its maximum, we have

and
or
so that in both cases

This is a contradiction to (8.35) and thus shows our claim.
By considering the functions, it follows, in particular, that a solution of the homogeneous heat equation attains both its maximum and its minimum on, so that in this case, with,
,
For the inhomogeneous equation one may show the following inequality, the proof of which we leave as an exercise, see Problem 8.7.
Theorem 8.7. The solution of (8.34) satisfies

Where r is the radius of a ball containing.
As usual such a result shows uniqueness and stability for the initial -boundary value problem.
We close this section by proving the uniqueness of a bounded solution to the pure initial value problem considered in Sect. 8.1.
Theorem 8.8. The initial value problem (8.1) has at most one solution which is bounded in, where T is arbitrary.
Proof. If there were two solutions of (8.1), then their difference would be a solution with initial data zero. It suffices therefore to show that the only bounded solution u of

is u = 0, or that and if is an arbitrary point in,and is arbitrary, then .We introduce the auxiliary function

and note that .Let now

Then

Since u is bounded we have on for some M. Defining R by, we have
,if,
And
,for
Hence we may apply Theorem 8.6 with and conclude that for .In particular, at we have, which completes the proof of the theorem.
The assumption of Theorem 8.8 that the solutions are bounded in may be relaxed to the requirement that for all, and for some M, c > 0, but without some such restriction on the growth of the solution for large , uniqueness is not guaranteed. For instance, the following function is a solution of the homogeneous heat equation which has initial v