实验１面向方程的数值积分方法仿真
1. 实验目的
运用CSS01.C 仿真程序解题,培养阅读及修改仿真程序的能力,学习并了解仿真程序的结构及特点.通过实验,加深理解4阶龙格-库塔法的原理及其稳定域.
2. 实验设备：装有BC 语言的PC 机一台
3. 实验内容
修改CSS01.C 仿真程序,对如下系统进行仿真.
(1) 线性定常系统
??????????321x x x =?????????
?---00600120000120??????????321x x x +??????????60000u y=[]00
1[]T
x x x 32
1
??????????)0()0()0(321x x x =????
?
?????000, u=1(t) a)、对龙格库塔法进行分析：它是一种数值积分法，也就是微分方程初
值问题数值计算法，是对初值微分方程的离散化求解。对于数值积分法我们常用的是欧拉法以及二阶和四阶龙格库塔法其原理分别如下： 对于形如
dt
dy
=),(y t f 的微分方程 用欧拉法仿真其迭代公式为 ),(*1k k k k y t f h y y +=+，其中h 为仿真步长（下同）；
二阶龙格库塔法其迭代公式为 y k+1= y k +h/2*(k1+k2)，其中),(1k k y t f k =，
);*,(21h k y h t f k k k ++=
四阶龙格库塔法其迭代公式为：)432221(*6/1k k k k h y y k k ++++=+，
其中),(1k k y t f k =，k2= f(t k +h/2,y k +k1*h/2)，k3= f(t k +h/2,y k +k2*h/2)，
k4= f(t k +h,y k +k3*h)；
本程序中大致分以下四次计算： 第一次计算：
g[1]*h/2——k[1][2]，y[1]——k[1][1]，把y[1]的值赋给k[1][1] Y[1]= k[1][1]+ k[1][2]；
dt
dy
=),(y t f ——g[1] ， （计算k1）， g[1]*h/2——k[1][2]，（把k1*h/2存储到k[1][2]） 第二次计算：
Y[1]= k[1][1]+ k[1][3]；
dt
dy
=),(y t f ——g[1] ， （计算k2）， g[1]*h/2——k[1][3]，（把k2*h/2存储到k[1][3]） 第三次计算：
dt
dy
=),(y t f ——g[1]，（计算k3）， Y[1]= k[1][1]+ k[1][4]； k[1][4]= g[1]*h ；（把k3*h 存储到k[1][4]） 第四次计算：
dt
dy
=),(y t f ——g[1]，（计算k4）， 此时 把k1,k2,k3,k4 的值代入下式：
)
432221(*6/1k k k k h y y k k ++++=+ 即为
6/)4*3**22**21*(1k h k h k h k h y y k k ++++=+ ，
由于h*k1=2*k[1][2]，2*h*k2=4*k[1][3]，2*h*k3=2*k[1][4]，h*k4=h*g[1] 所以可得如程序中的计算公式：
Y[1]= k[1][1]+(2*k[1][2]+4* k[1][3]+2* k[1][4]+h*g[1]）； 可见该仿真程序为四阶龙格库塔法的数值积分仿真程序。
以上为只是一个微分方程的情况，若是微分方程组其过程也大概如此。只是将数组y 、g 和k 的其它部分也相应赋值或做运算即可。其实是把g[1],y[1]分别用g[i],y[i]替代；k[1][j]用k[i][j]替代即可（j=1、2、3、4）。 b)、计算)(t y 理论表达式，首先引入状态变量: y x =1
dt dy x
x ==12 2223dt y d x x == 则有
y dt dy x 202+=
222320020x dt dy dt y d x ++=
u y dt dy dt y d dt y d 60060040002202233=+++
600
4000220600
)(23+++=
∴s s s s G )
6125.61)(05.01)(005.01(1
)15125.0)(9832.19)(0165
.200(600
600400220600)(23s s s s s s s s s s +++=
+++≈
+++=
φ 则)
6125.61)(05.01)(005.01(1
)(s s s s s Y +++=利用分式部分法，取拉氏
反变换有 s
A s A s A s A 6125.6105.01005.014
321++++++=
1)(01=?==s s s Y A 7165.200210203.4)005.01()(-=?-=+?=s s Y A 49832.1931024.4)05.01()(--=?=+?=s s s Y A 62158.6)6125.61)((15125.04-=+=-=s s s Y A
t
t t e e e t y 0165.20079832.19415125.010204.41024.462158.61)(-----?-?+-=∴ 0≥t
c)、程序实现：在CSS01.C 仿真程序基础上,增加float u ； i n t 1=u 将
]6][102[w 再修改为]6][801[w 在
100
323(}100
33({
><> 中的将100改 为800，再在);101 ;1(++<==i i i for 中将101改为801 最后在输入程序块)(mod sub 程序改为； ; *600]1[*600]3[];3[]2[*200]2[];2[]1[*20]1[]; 1[]22[;]21[u y g y y g y y g y y u y +-=+-=+-=== 在输出转换)(repsub 程序块中程序改为； ];21[]1[y a = ];22[]2[y a =];2[]3[y a = ];3[]4[y a = 以下确定仿真步长：由上面方程可得系统传递函数为 600 )20)(200(600 )(+++= s s s s G 其开环惯量为201+s ， 2001+s 理论上仿真步长h<2|1| a (a 为惯量极点)，则有h<2|200 1 |=0.01，可取h 为0.01。由于程序中存储的仿真点数最多为800个，所以可取 总时间为8。 通过作以上的修改，等程序运行后分别输入T1=6.6,T2=0.01 N1=3 N2=0 N3=661 J8=1 Input initial values of state variables=0 回车 0回车 0回车 实验结果：由仿真结果可以看出：该模型是一个系统阶跃响应实 例。 当t 趋于无穷时y(t)的理论值趋于1若取t 为k*0.01(k=1,2,3,4,5,6…)计算y(t)的值并与仿真结果中的y[1]进行比较，发现其值大致相等，所以该仿真是正确的，合理的。(程序见附录xianxin) 实验小结：再用根轨迹法分析： 1）、确定实轴上的根轨迹，实轴上[0，—200.0165]区域必为根轨迹 2）、确定渐进线（根轨迹） 3=-m n 故有三条轨迹渐进线 333.73-=a σ 60)12(=-+= m n k a π ? ， 180 ， 300 （ k 分 别为0、1、２ ） ３）、求分离点：由 01 1=-∑=n i i p d 有 0200 1 20101=++++-d d d 既转化为 0400044032=++d d 7373865722 .9112-==∴d s 92928.1363-=s （舍去）设分离点处开环增益为k 把 7373865722.9112-==d s 代入04000220 2 3=+++k s s s 中得 13969.19013=k 当 13969.19013=k 时临界阻尼 当 13969.19013>k 时为欠阻尼 4）、确定根轨迹与虚轴交点：本题闭环特征方程式为： 04000220*3=+++k s s s 对上式应用劳斯判据，有 * * 1*2 30 2204000220220 40001k s k s k s s -? 令劳斯表中1s 行首项为0 得880000* =k 根据2 s 行的系数，得如下辅助方程：0220* 2 =+k s 代入880000* =k 并令jw s =解出交点坐标2455532.63±=w 。 所以当超过这两点系统进入不稳定状态。 （2） 实验内容----非线性系统（弱肉强食模型）： ???????+-=-=)()()()()()()() (t y t bx t sy dt t dy t y t ax t rx dt t dx (1) 其中:001.0=r ，6 102-?=a ，12000)0(=x ，01.0=s ， 6101-?=b ，600)0(=y 。 a)、 理论分析： 设兔子数量为)(t x ,老虎的数量为)(t y ,将它们置于一个 特定的环境里, )(t x 将以自然增长率r 增长,即rx x = .但是老虎以兔子为食,致使兔子的增长率降低,设降低程度与老虎数量)(t y 成正比,于是相对增长率为ay r r x -=.常数a 反映了老虎捕食兔子的能力.如果没有兔子,老虎就 无法生存,设老虎的自然死亡率为s ,则sy y -= .兔子为老虎提供了食物,致使老虎死亡率降低,或者说为老虎提供了增长的条件.设增长率与兔子的数量 )(t x 成正比,于是老虎的相对增长率bx s r y +-=.常数b 反映了兔子对老虎 的供养能力.所以形成上面的模型. 下面对上述非线性微分方程组做稳定性分析. 首先求平衡点,令,0,0==y x 得(0,0)和),(a r b s .若只考虑),(a r b s 的稳定性,则可作变量代换，使原点为平稳点,令) ()(),()(t a r t y t b s t x ψρ+=+=代入(1)式得???????+=+-=) ()(ψρψρψρ a r b dt d b s a dt d 设b s <<ρ,a r <<ψ,去掉二次项使之线 性化,得???????=-=ρψψρa br dt d b as dt d 即?????????? ? ? ????- =??????????ψρψρ00a br b as dt d dt d (2) 可以求出上式系数行列式的特征根为rsi ±.故线性方程组的通解为 ?? ? ??==)sin()()cos()(00rsi t rsi t ψψρρ (3) 由上式可以看出: (a)ψρ,的线性化解,既不增长也不衰减,而是连续振动.这意味着平稳是亚稳定的,这是一种广义稳定(Kolmogorov 意义下);在平衡点领域显示出稳定的有界循环. (b) 令π2=rsT ,则振动周期T=rs π 2 (c) 兔子和老虎的总数的振动相位差为0 90 ,见 图 兔子和老虎数量变化图示 (4) 从解中消去时间t,得12 02 =??? ? ??+???? ??ψψρ ρ 这 条 轨 线 在 总 体 相 空 间 的 图 形 如 图 轨线在总体相空间的图形 （横虚线表示λr ） (5) 可定出在相空间运动的方向,将)()(),()(t a r t y t b s t x ψρ+=+= 代入线性化方程可得??? ? ??? -=--=)()(b s x a br dt dy a r y b as dt dx 因此轨线在总体相空间的图形四个区域中有 I:;0,0> II:;0,0< x(t)和 y(t)的平均值,则 b s dt y y s b T dt t x T x T =+==??)(11)(10 其中y(t)为周期函数,y(T)=y(0). 同理,a r y = . 对于周期性变化的x(t)和y(t),用它们的平均值y x ,表示其大小.上面的 分析表明,兔子的平均数量取决于老虎方程)(bx s y y +-= 中的参数s 和b,而老虎的平均数量取决于兔子方程)(ay r x x -= 中的参数r 和a ．当兔子的 自然增长率r 下降时，老虎的数量将减少，而当老虎捕食兔子的能力提高时，对兔子没有影响，只使老虎减少．另一方面，老虎死亡率s 上升将导致兔子增多，而兔子对老虎的供养能力b 的提高，会导致兔子减少。 b)、程序实现：同样在CSS01.C 仿真程序基础上作如下修改，先将 100 323(}100 33({ ><> ;101;1(++<==i i i for 中将101改为801最后在系统模型输入程序块)(mod sub 中程序改为； ]; 2[*]1[*000001.0]2[*01.0]2[]; 2[*]1[*000002.0]1[*001.0]1[y y y g y y y g +-=-= 在输入转换)(repsub 程序块中程序改为； ];1[]1[y a = ];2[]2[y a =];1[]3[g a = ];2[]4[g a = 因为在这个封闭的生态系统里，循环周期较长，故可取仿真总时间为６４００天，步长为１０天，即T1=6400,T2=10 N1=2 N2=0 N3=640 J8=1 Input initial values of state variables=12000回车 600回车 0回车 作图程序： void ()blot { double ;1;0==j i );"\\\\:",mod ,&(&bgi bc d e g gdrive initgraph );2000 (delay );800 ;0(++<==i i i for { );(RED setcolor ,(i i moveto ); (int)( W [ /]3 50 ][ (int)(+ i lineto [ W /]3 1 ) 50 ][ ); setcolor (WHITE ); ,(i W moveto+ 200 i [ (int)( 5/]2 ][ )); ,(+ W lineto 200 i +i ); (int)( 1 )5/]2 [ ][ } getch (); d bc (&bgi g gdrive initgraph e " ); ,& \\ \\: mod ", ( delay 2000 ); pr(" "x int y change with); f t h e s =s fpr <= s ) ; 800 ;0 (+ + { (int)(- [ =s W i 10 ][ ; 700 ) /]3 W j= (int)(s 3/]2 ); ][ [ (GREEN setcolor ); moveto i ); ,(j i lineto j ,(+ 1 ); } getch (); } 实验结果：仿真结果与理论分析图形基本一样如前面两图所示：兔 子和老虎的数量各自在趋于一个动态平衡。即在一定的值上下来回振动，这是一种广义的平衡状态。它们之间数量关系近似为一个椭圆 （程序见附录feixian1） 4、实验小结：对数值积分的稳定性分析：一个本来 是稳定的系统，利用数值积分法进行仿真的时候常常会得出不稳定的结论。造成这种现象的原因往往是计算的字长选取的太大，当步h 选取的太大的时候，数值积分法会使得各种误差传递出去，从而引起计算不稳定。而误差的来源一般是以下几个：初始误差（在实际计算中给出的初值0y 通常会不太准确）；由于计算机字长限制造成的舍入误差：在某一步长h 下产生的截断误差。而这些误差都可能在计算中向下传播。 数值解的稳定性，是指在拢动（初始误差、舍入、截断误差等效影响下，其计算过程中的累积误差不会随计算步数的增加而无限增长。不同的数值法对应着不同的差分递推公式。一个数值是否稳定取决于该差分方程的特征根是否满足稳定性要求。 下面以一阶显式亚当姆斯方法来讨论数值积分的稳定性分析方法 假定系统的微分方程为 y dt dy λ=， βαλi += ，0<α （１） 其递推公式为 n n n y h y y λ+=+1 （２） 假定k y （k=0,1,2,…..）为它的一个仿真解，另 外 设 ) (k k y ε+为其准确解，则有 )()()(11n n n n n n y h y y ελεε+++=+++ （３） （３）－（２）得 n n n h λεεε+=+1 整理后为 0)1(1=+-+n n h ελε 其特征方程为0)1(=+-λh z 为了使误差 序列k ε不随k 的增加而增加，必须要求其特征根λh z +=1在单位圆内，即11≤+λh ，它所对应的区域就是一阶显式亚当姆斯算法的稳定域。 因此，如果系统方程0,<=ααλ，运用一阶显式亚当姆斯算法时，则只有当α12≤h 才能保证计算是稳定的。 电子一班王申江 实验九数值微积分与方程数值求解 一、实验目的 1、掌握求数值导数和数值积分的方法 2、掌握代数方程数值求解的方法 3、掌握常微分方程数值求解的方法 二、实验内容 1、求函数在指定点的数值导数。 () 23 2 123,1,2,3 026 x x x f x x x x x == >>syms x >>f=[x x^2 x^3;1 2*x 3*x^2;0 2 6*x]; >>F=det(f) F=2*x^3 >>h=0.1 >>x=[0:h:4]; >>f=2*x^3; >>[dy,dx]=diff_ctr(f,h,1); >>y1=dy(dx==1) y1=6.0000 >>y2=dy(dx==2) y2=24.0000 >>y3=dy(dx==3) y3=54.0000 2、用数值方法求定积分。 （1） 210I π =?的近似值 a=inline('sqrt(cos(t.^2)+4*sin((2*t).^2)+1)'); I=quadl(a,0,2*pi) I = 6.7992 + 3.1526i （2）()1 202ln 11x I dx x +=+? b=inline('log(1+x)./(1+x.^2)'); I=quadl(b,0,1) I = 0.2722 3、分别用3种不同的数值方法解线性方程组。 6525494133422139211 x y z u x y z u x y z u x y u +-+=-??-+-=??++-=??-+=? A=[6,5,-2,5;9,-1,4,-1;3,4,2,-2;3,-9,0,2]; b=[-4,13,1,11]'; x=A\b 几种定积分的数值计算方法 摘要:本文归纳了定积分近似计算中的几种常用方法,并着重分析了各种数值方法的计 算思想,结合实例,对其优劣性作了简要说明. 关键词:数值方法;矩形法;梯形法;抛物线法;类矩形;类梯形 Several Numerical Methods for Solving Definite Integrals Abstract:Several common methods for solving definite integrals are summarized in this paper. Meantime, the idea for each method is emphatically analyzed. Afterwards, a numerical example is illustrated to show that the advantages and disadvantages of these methods. Keywords:Numerical methods, Rectangle method, Trapezoidal method, Parabolic method, Class rectangle, Class trapezoid 1. 引言 在科学研究和实际生产中,经常遇到求积分的计算问题,由积分学知识可知,若函数 )(x f 在区间],[b a 连续且原函数为)(x F ,则可用牛顿-莱布尼茨公式 ?-=b a a F b F x f ) ()()( 求得积分.这个公式不论在理论上还是在解决实际问题中都起到了很大的作用. 在科学研究和实际生产中,经常遇到求积分的计算问题,由积分学知识可知,若函数)(x f 在区间],[b a 连续且原函数为)(x F ,则可用牛顿-莱布尼茨公式 ?-=b a a F b F x f ) ()()( 求得积分.这个公式不论在理论上还是在解决实际问题中都起到了很大的作用.另外,对于求导数也有一系列的求导公式和求导法则.但是,在实际问题中遇到求积分的计算,经常会有这样的情况: (1)函数)(x f 的原函数无法用初等函数给出.例如积分 dx e x ?-1 02 , ? 1 sin dx x x 等,从而无法用牛顿-莱布尼茨公式计算出积分。 (2)函数)(x f 使用表格形式或图形给出,因而无法直接用积分公式或导数公式。 (3)函数)(x f 的原函数或导数值虽然能够求出,但形式过于复杂,不便使用. 由此可见,利用原函数求积分或利用求导法则求导数有它的局限性,所以就有了求解数值积分的很多方法,目前有牛顿—柯特斯公式法,矩形法,梯形法,抛物线法,随机投点法,平均值法,高斯型求积法,龙贝格积分法,李查逊外推算法等等,本文对其中部分方法作一个比较． 2.几何意义上的数值算法 s 在几何上表示以],[b a 为底,以曲线)(x f y =为曲边的曲边梯形的面积A ,因此,计 算s 的近似值也就是A 的近似值,如图1所示.沿着积分区间],[b a ,可以把大的曲边梯形分割成许多小的曲边梯形面积之和.常采用均匀分割,假设],[b a 上等分n 的小区间 ,x 1-i h x i +=b x a x n ==,0,其中n a b h -= 表示小区间的长度. 2.1矩形法 实验二数值方法计算积分 学号：姓名：指导教师：实验目的 1、了解并掌握matlab软件的基本编程、操作方法； 2、初步了解matlab中的部分函数，熟悉循环语句的使用； 3、通过上机进一步领悟用复合梯形、复合辛普森公式，以及用龙贝格求积 方法计算积分的原理。 一、用不同数值方法计算积分 10x ln xdx=-94. (1)取不同的步长h.分别用复合梯形及辛普森求积计算积分，给出误差中关 于h的函数，并与积分精确值比较两个公式的精度，是否存在一个最小 的h，使得精度不能再被改善？ (2)用龙贝格求积计算完成问题（1）。 二、实现实验 1、流程图： 下图是龙贝格算法框图： 2、 算法： (1) 复合梯形公式：Tn=++)()([2b f a f h 2∑-=1 1 )](n k xk f ； (2) 复合辛普森公式：Sn=6h [f(a)+f(b)+2∑-=11)](n k xk f +4∑-=+1 )2/1(n k x f ]； 以上两种算法都是将a-b 之间分成多个小区间（n ）,则h=(b-a)/n,x k =a+kh, x k+1/2=a+(k+1/2)h,利用梯形求积根据两公式便可。 (3) 龙贝格算法：在指定区间内将步长依次二分的过程中运用如下公式 1、Sn= 34T2n-31 Tn 2、 Cn=1516S2n-151 Sn 3、 Rn=6364C2n-631 Cn 从而实现算法。 3、 程序设计 （1）、复合梯形法： function t=natrapz(fname,a,b,n) h=(b-a)/n; fa=feval(fname,a);fb=feval(fname,b);f=feval(fname,a+h:h:b-h+0. 001*h); t=h*(0.5*(fa+fb)+sum(f)); （2）、复合辛普森法： function t=natrapz(fname,a,b,n) h=(b-a)/n; fa=feval(fname,a);fb=feval(fname,b);f1=feval(fname,a+h:h:b-h+0 .001*h); f2=feval(fname,a+h/2:h:b-h+0.001*h); t=h/6*(fa+fb+2*sum(f1)+4*sum(f2)); （3）龙贝格法： function [I,step]=Roberg(f,a,b,eps) if(nargin==3) eps=1.0e-4; end; M=1; tol=10; k=0; T=zeros(1,1); h=b-a; T(1,1)=(h/2)*(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)), 《微分方程数值解法》 【摘要】自然界与工程技术中的很多现象，可以归结为微分方程定解问题。其中，常微分方程求解是微分方程的重要基础内容。但是，对于许多的微分方程，往往很难得到甚至不存在精确的解析表达式，这时候，数值解提供了一个很好的解决思路。，针对于此，本文对常微分方程数值解法进行了简单研究，主要讨论了一些常用的数值解法，如欧拉法、改进的欧拉法、Runge —Kutta 方法、Adams 预估校正法以及勒让德谱方法等，通过具体的算例，结合MA TLAB 求解画图，初步给出了一般常微分方程数值解法的求解过程。同时，通过对各种方法的误差分析，让大家对各种方法的特点和适用范围有一个直观的感受。 【关键词】 常微分方程 数值解法 MA TLAB 误差分析 引言 在我国高校，《微分方程数值解法》作为对数学基础知识要求较高且应用非常广泛的一门课程，不仅 在数学专业，其他的理工科专业的本科及研究生教育中开设这门课程．近四十年来，《微分方程数值解法》不论在理论上还是在方法上都获得了很大的发展．同时，由于微分方程是描述物理、化学和生物现象的数学模型基础，且它的一些最新应用已经扩展到经济、金融预测、图像处理及其他领域 在实际应用中，通过相应的微分方程模型解决具体问题，采用数值方法求得方程的近似解，使具体问题迎刃而解。 2 欧拉法和改进的欧拉法 2.1 欧拉法 2.1.1 欧拉法介绍 首先，我们考虑如下的一阶常微分方程初值问题 ???==0 0)() ,('y x y y x f y （2--1） 事实上，对于更复杂的常微分方程组或者高阶常微分方程，只需要将x 看做向量，（2--1）就成了一个一阶常微分方程组，而高阶常微分方程也可以通过降阶化成一个一阶常微分方程组。 欧拉方法是解常微分方程初值问题最简单最古老的一种数值方法，其基本思路就是把（2--1）中的导数项'y 用差商逼近，从而将一个微分方程转化为一个代数方程，以便求解。 设在[]b a ,中取等距节点h ，因为在节点n x 点上，由（2--1）可得： 1. 课本2p 有证明 2. 课本812,p p 有说明 3. 课本1520,p p 有说明 4. Rit2法，设n u 是u 的n 维子空间，12,...n ???是n u 的一组基底，n u 中的任一元素n u 可 表为1n n i i i u c ?==∑ ，则,11 11()(,)(,)(,)(,)22j n n n n n n i j i j j i j j J u a u u f u a c c c f ???=== -=-∑∑是12,...n c c c 的二次函数，(,)(,)i j j i a a ????=，令 () 0n j J u c ?=?，从而得到12,...n c c c 满足1 (,)(,),1,2...n i j i j i a c f j n ???===∑，通过解线性方程组，求的i c ，代入1 n n i i i u c ?==∑， 从而得到近似解n u 的过程称为Rit2法 简而言之，Rit2法：为得到偏微分方程的有穷维解，构造了一个近似解，1 n n i i i u c ?== ∑， 利用,11 11()(,)(,)(,)(,)22j n n n n n n i j i j j i j j J u a u u f u a c c c f ???===-=-∑∑确定i c ，求得近似解n u 的过程 Galerkin 法：为求得1 n n i i i u c ? == ∑形式的近似解，在系数i c 使n u 关于n V u ∈，满足(,)(,) n a u V f V =，对任 意 n V u ∈或（取 ,1j V j n ?=≤≤） 1 (,)(,),1,2...n i j i j i a c f j n ???===∑的情况下确定i c ，从而得到近似解1 n n i i i u c ?==∑的过程称 Galerkin 法为 Rit2-Galerkin 法方程： 1 (,)(,)n i j i j i a c f ???==∑ 5. 有限元法：将偏微分方程转化为变分形式，选定单元的形状，对求解域作剖分，进而构 造基函数或单元形状函数，形成有限元空间，将偏微分方程转化成了有限元方程，利用 有效的有限元方程的解法，给出偏微分方程近似解的过程称为有限元法。 6. 解：对求解区间进行网格剖分，节点01......i n a x x x x b =<<<<=得到相邻节点1,i i x x - 第十五章 积分方程 积分方程论是泛函分析的一个重要分支，它是研究数学其他学科（例如偏微分方程边值问题）和各种物理问题的一个重要数学工具。本章叙述线性积分方程，重点介绍弗雷德霍姆积分方程的性质和解法；并简略地介绍了沃尔泰拉积分方程以及一些奇异积分方程；此外，还扼要地叙述积分方程的逐次逼近法和预解核，并举例说明近似解法；最后考察了一个非线性积分方程。 §1 积分方程一般概念与弗雷德霍姆方程 一. 积分方程一般概念 1. 积分方程的定义与分类 [线形积分方程] 在积分号下包含未知函数y (x )的方程 ()()()()(),d b a x y x F x K x y αλξξξ=+? (1) 称为积分方程。式中α(x ),F (x )和K (x,ξ)是已知函数，λ,a,b 是常数，变量x 和ξ可取区间(a,b ) 内的一切值；K (x,ξ)称为积分方程的核，F (x )称为自由项，λ称为方程的参数。如果K (x,ξ)关于x,ξ是对称函数，就称方程(1)是具有对称核的积分方程；如果方程中的未知函数是一次的，就称为线性积分方程，方程(1)就是线性积分方程的一般形式；如果F (x )≡0 ，就称方程(1)为齐次积分方程，否则称为非齐次积分方程。 [一维弗雷德霍姆积分方程（Fr 方程）] 第一类Fr 方程 ()()(),d b a K x y F x ξξξ=? 第二类Fr 方程 ()()()(),d b a y x F x K x y λξξξ=+? 第三类Fr 方程 ()()()()(),d b a x y x F x K x y αλξξξ=+? [n 维弗雷德霍姆积分方程] 111()()()()(),d D P y P F P K P P y P P α=+? 称为n 维弗雷德霍姆积分方程，式中D 是n 维空间中的区域，P ,P 1∈D ，它们的坐标分别是 (x 1,x 2, ,x n )和),,,(21 n x x x ''' ,α(P )=α(x 1,x 2, ,x n ),F (P )=F (x 1,x 2, x n )和K (P ,P 1)=K (x 1,x 2, ,x n , ),,,21 n x x x ''' 是已知函数，f (P )是未知函数。 关于Fr 方程的解法，一维和n (>1)维的情况完全类似，因此在以后的讨论中仅着重考虑一维Fr 方程。 [沃尔泰拉积分方程] 如果积分上限b 改成变动上限，上面三类Fr 方程分别称为第一、第二、第三类沃尔泰拉积分方程。 由于第三类Fr 方程当α(x )在(a ,b )内是正函数时，可以化成 求第一类Fredholm积分方程的离散正则化方法 【摘要】基于矩阵奇异值分解的离散正则化算法，本文给出了第一类Fredholm积分方程的求解方法。并通过算例验证了此算法的可行性。 【关键词】第一类Fredholm积分方程；矩阵奇异值；正则化方法 0 引言 在实际问题中，有很多数学物理方程反问题的求解最后总要归结为一个第一类算子方程： Kx=y（1） 的求解问题，其中K是从Hilbert空间X到Hilbert空间Y一个有界线性算子，x∈X，y∈Y。通常右端项y是观测数据，因而不可避免的带有一定的误差δ。文中假设方程（1）的右端的扰动数据yδ∈Y满足条件：yδ-y≤δ（C1）。我们需要求解扰动方程Kx=yδ∈Y。（2） 通常境况下，当K为紧算子时，方程（1）的求解时不适定的[1]。即右端数据的小扰动可导致解的巨大变化。消除不稳定性的一个自然的方式是用一族接近适定问题的模型去逼近原问题，比如最著名的Tikhonov正则化方法，用如下适定的算子方程： 去逼近原问题Kx=yδ，其中α>0为一正的“正则参数”，K*表示K的伴随算子。正则化[2-3]是近似求解方程（1）的一种有效方法。Krish应用奇异系统理论提出的正则化子的概念，这给正则化方法的建立提供了新的理论依据。本文利用基于矩阵奇异值分解的离散正则化算法，通过适当选取正则化参数进行不适定问题的求解。 1 基于矩阵奇异值分解的离散正则化算法 矩阵的奇异值分解（SVD）是现代数值线性代数中最重要的基本计算分析工具之一，它具有优良的数值稳定性。其重要应用领域包括矩阵理论以及自动控制理论，力学和物理学等，还有更多的应用方面尚在继续探索中。 对于一般算子方程Kx=y，利用高斯-勒让德求积公式、复化梯形公式或者复化辛普森求积公式等的数值方法将它离散得到一个矩阵方程Ax=y，这样，算子方程Kx=y的求解就转化为矩阵方程： 的求解。 定义设A是m×n实矩阵（m≥n），称n阶方阵ATA的非零特征值的算术平 目录 第一章数值积分计算的重述 (1) 1.1引言 (1) 1.2问题重述 (2) 第二章复化梯形公式 (3) 2.1 复化梯形公式的算法描述 (3) 2.2 复化梯形公式在C语言中的实现 (3) 2.3 测试结果 (4) 第三章复化simpson公式 (6) 3.1 复化simpson公式的算法描述 (6) 3.2 复化simpson公式在C语言中的实现 (6) 3.3 测试结果 (7) 第四章复化cotes公式 (8) 4.1 复化cotes公式的算法描述 (8) 4.2 复化cotes公式在C语言中的实现 (9) 4.3 测试结果 (10) 第五章Romberg积分法 (11) 5.1 Romberg积分法的算法描述 (11) 5.2 Romberg积分法在C中的实现 (12) 5.3 测试结果 (13) 第六章结果对比分析和体会 (144) 参考文献 (16) 附录 (16) 数值积分?-10 2 dx e x （一） 第一章 数值积分计算的重述 1.1引言 数值积分是积分计算的重要方法，是数值逼近的重要内容，是函数插值的最直接应用，也是工程技术计算中常常遇到的一个问题。在应用上，人们常要求算出具体数值，因此数值积分就成了数值分析的一个重要内容。在更为复杂的计算问题中，数值积分也常常是一个基本组成部分。 在微积分理论中，我们知道了牛顿-莱布尼茨(Newton-Leibniz)公式 ()() () b a f x d x F b F a =-? 其中()F x 是被积函数()f x 的某个原函数。但是随着学习的深入，我们发现一个问题: 对很多实际问题，上述公式却无能为力。这主要是因为:它们或是被积函数没有解析形式的原函数,或是只知道被积函数在一些点上的值,而不知道函数的形式,对此,牛顿—莱布尼茨(Newton-Leibniz)公式就无能为力了。此外,即使被积函数存在原函数,但因找原函数很复杂,人们也不愿花费太多的时间在求原函数上,这些都促使人们寻找定积分近似计算方法的研究,特别是有了计算机后,人们希望这种定积分近似计算方法能在计算机上实现,并保证计算结果的精度,具有这种特性的定积分近似计算方法称为数值积分。由定积分知识,定积分只与被积函数和积分区间有关,而在对被积函数做插值逼近时,多项式的次数越高,对被积函数的光滑程度要求也越高,且会出现Runge 现象。如7n >时,Newton-Cotes 公式就是不稳定的。因而,人们把目标转向积分区间,类似分段插值,把积分区间分割成若干小区间,在每个小区间上使用次数较低的Newton-Cotes 公式,然后把每个小区间上的结果加起来作为函数在整个区间上积分的近似,这就是复化的基本思想。本文主要 微分方程的分类及其数值解法 微分方程的分类： 含有未知函数的导数，如dy/dx=2x 、ds/dt=0.4都是微分方程。 一般的凡是表示未知函数、未知函数的导数与自变量之间的关系的方程，叫做微分方程。未知函数是一元函数的，叫常微分方程；未知函数是多元函数的叫做偏微分方程。微分方程有时也简称方程。 一、常微分方程的数值解法： 1、Euler 法： 00d (,), (1.1)d (), (1.2) y f x y x y x y ?=???=? 001 (),(,),0,1,,1n n n n y y x y y hf x y n N +=??=+=-? (1.4) 其中0,n b a x x nh h N -=+=. 用(1.4)求解(1.1)的方法称为Euler 方法。 后退Euler 公式???+==+++),,(),(111 00n n n n y x hf y y x y y 梯形方法公式 )].,(),([2 111+++++=n n n n n n y x f y x f h y y 改进的Euler 方法11(,),(,),1().2p n n n c n n p n p c y y hf x y y y hf x y y y y ++?=+??=+???=+??? 2、Runge-Kutta 方法： p 阶方法 ： 1()O h -=?总体截断误差局部截断误差 二阶Runge-Kutta 方法 ??? ????++==++=+),,(),,(,2212 1211hk y h x f k y x f k k h k h y y n n n n n n i.常微分方程初值问题数值解法 常微分方程初值问题的真解可以看成是从给定初始点出发的一条连续曲线。差分法是常微分方程初值问题的主要数值解法，其目的是得到若干个离散点来逼近这条解曲线。有两个基本途径。一个是用离散点上的差商近似替代微商。另一个是先对微分方程积分得到积分方程，再利用离散点作数值积分。 i.1 常微分方程差分法 考虑常微分方程初值问题：求函数()u t 满足 (,), 0du f t u t T dt =<≤ （i.1a ） 0(0)u u = (i.1b) 其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的连续函数，0u 和T 是给定的常数。我们假设(,)f t u 对u 满足Lipschitz 条件，即存在常数L 使得 121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-?∈∈-∞∞ （i.2） 这一条件保证了（i.1）的解是适定的，即存在，唯一，而且连续依赖于初值0u 。 通常情况下，（i.1）的精确解不可能用简单的解析表达式给出，只能求近似解。本章讨论常微分方程最常用的近似数值解法-差分方法。先来讨论最简单的Euler 法。为此，首先将求解区域[0,]T 离散化为若干个离散点： 0110N N t t t t T -=<< <<= （i.3） 其中n t hn =，0h >称为步长。 在微积分课程中我们熟知，微商（即导数）是差商的极限。反过来，差商就是微商的近似。在0t t =处，在（i.1a ）中用向前差商 10()()u t u t h -代替微商du dt ，便得 10000()()(,())u t u t hf t u t ε=++ 如果忽略误差项0ε，再换个记号，用i u 代替()i u t 便得到 1000(,)u u hf t u -= 一般地，我们有 1Euler (,), 0,1, ,1n n n n u u hf t u n N +=+=-方法： （i.4） 从(i.1b) 给出的初始值0u 出发，由上式可以依次算出1,,N t t 上的差分解1,,N u u 。 一类第二种非线性Volterra 积分方程 积分数值解方法 1前言 微分方程和积分方程都是描述物理问题的重要数学工具，各有优点.相对于某种情况来说，对于某种物理数学问题，积分方程对于问题的解决比微分方程更加有优势，使对问题的研究更加趋于简单化，在数学上，利用积分形式讨论存在性、唯一性往往比较方便，结果也比较完美,所以研究积分方程便得越来越有用,日益受到重视. 积分方程的发展，始终是与数学物理问题的研究息息相关.一般认为，从积分发展的源头可以追溯到国外的数学家克莱茵的著作《古今数学思想》，该书是被认为第一个清醒的认为应用积分方程求解的是Abel.Abel 分别于1833年和1826年发表了两篇有关积分方程的文章，但其正式的名称却是由数学家du Bois-Raymond 首次提出的，把该问题的研究正式命名为积分方程。所以最早研究积分方程的是Abel,他在1823年从力学问题时首先引出了积分方程，并用两种方法求出了它的解,第一的积分方程便是以Abel 命名的方程.该方程的形式为：?=-b a a x f dt t x t )()() (?,该方程称为广义Abel 方程,式中a 的值在(0,1)之间.当a=21时，该式子便成为)()(x f dt t x x x a =-??.在此之前，Laplace 于1782年所提出的求Laplace 反变换问题，当时这个问题就要求解一个积分方程.但是Fourier 其实已经求出了一类积分方程的反变换，这就说明在早些时候积分方程就已经在专业性很针对的情况下得到了研究，实际上也说明了Fourier 在研究反变换问题是就相当于解出了一类积分方程.积分方程的形成基础是有两位数学家Fredholm 和V olterra 奠定的，积分方程主要是研究两类相关的方程，由于这两位数学家的突出贡献，所以这两个方程被命名为Fredholm 方程和V olterra 方程。后来又有德国数学家D.Hilbert 进行了重要的研究，并作出了突出的贡献，由于D.Hilbert 领头科学家的研究，所以掀起了一阵研究积分方程的热潮，并出现了很多重要的成果，后来该理论又推广到非线性部分。我国在60年代前，积分方程这部分的理论介绍和相关书本主要靠翻译苏联的相关书籍，那时研究的积分方程基本是一种模式，即用古典的方法来研究相关的积分方程问题，这样使得问题的研究变得繁琐、复杂，在内容方面比较单一、狭隘，甚至有些理论故意把积分方程的研究趋向于复杂化。随着数学研究的高速发展，特别是积分方程近年来的丰富发展，如此单一、刻板的解法已经不能跟上数学研究时代的步伐。在九十年代我国的数学专家路见可、钟寿国出版了《积分方程论》，该书选择2L 空间来讨论古典积分方程，并结合泛函分析的算子理论来分析积分方程的相关问题。最近出版的比较适 第一篇积分方程 第一章方程的导出和基本概念 §1.1 方程的导出 许多力学、工程技术和数学物理问题都能用积分方程形式描述，而求解常微分方程和偏微分方程的定解问题常常可转化为求解积分方程的问题。下面举几个典型的问题作为例子，扼要地阐明导出积分方程的方法以及微分方程与积分方程之间的联系。 例1：弹性弦负荷问题 一根轻且软的弹性弦,长为l,两端固定,如图所示,静止时与x轴重合,弦内张力为 T.今在其上加以强度为 ()x ?的负荷.设在任一点M (横坐标 为x ) ()x ?, 且设 解：在任一点x ξ=处取微小的一段弦d ξ,则作用于其上的重力为 ()d ?ξξ,记之为0P ,则这一重力0P 必 引起弦的形变,记ξ处位移为S ,则： 01020sin sin T T P θθ+=, 因为0()T x ?>>,所以12,1θθ<< 112sin tan ,sin .S S l θθθξξ ?≈=≈- 所以000S S T T P l ξξ ?+? =-, 得 00()P l S T l ξξ-=?. 记0P 引起的x 处位移为* ()y x , 则0x ξ≤≤时, 由y S x ξ *=得 * 00() ()P l S y x x x T l ξξ-=?=??; 当x l ξ≤≤时，y S l x l ξ*= -- , ? 00()()P l x y x T l ξ* -= ??; 记：0 0,0(,),.l x x T l G x l x x l T l ξ ξξξξ-??≤≤??=?-??≤≤?? 则 0()(,)y x G x P ξ* =, ()(,)()y x G x d ξ?ξξ* =， 对ξ从0l 到求积分, 第七章 数值积分 如果函数f(x)在区间[a,b]上连续，且原函数为F(x)，则可用牛顿―莱布尼兹 公式：)()()(a F b F dx x f b a -=?来求得定积分。然而很多函数无法用牛顿―莱布尼兹公式求积分。 一个简单被积函数，例如，其不定积分可能很 复杂，见下面的MA TLAB 实例： >> syms a b c x >> int(sqrt(a+b*x+c*x*x),x) ans=1/4*(2*c*x+b)/c*(a+b*x+c*x^2)^(1/2)+1/2/c^(1/2)*log((1/2*b+c*x )/c^(1/2)+(a+b*x+c*x^2)^(1/2))*a-1/8/c^(3/2)*log((1/2*b+c*x)/c^(1/2)+(a+b*x+c*x^2)^(1/2))*b^2 所以有必要研究简单、高效的计算定积分的方法（即数值积分方法）。数值积分的基本思想是构造一个简单函数P n (x )来近似代替被积分函数f (x )，然后通过求?b a n dx x P )(得?b a dx x f )(的近似值。 7.1 插值型求积公式 设?=b a dx x f I )(*，插值型求积公式就是构造插值多项式P n (x )，使 ?=≈b a n dx x P I I )(* 。 构造以a ，b 为结点的线性插值多项式)()()(1b f a b a x a f b a b x x P --+--= ，[])()()(21)()()(1b f a f a b dx b f a b a x a f b a b x dx x P T b a b a +-=?? ? ???--+--==??称为梯形公式。 积分方程理论的发展，始终与数学物理问题的研究紧密相联，它在工程、力学等方面有着极其广泛的应用。通常认为，最早自觉应用积分方程并求出解的是阿贝尔(Abel)，他在1823年研究质点力学问题时引出阿贝尔方程。此前，拉普拉斯(Laplace)於1782年在数学物理中研究拉普拉斯变换的逆变换以及傅里叶(Fourier)於1811年研究傅里叶变换的反演问题实际上都是解第一类积分方程。随着计算技术的发展，作为工程计算的重要基础之一，积分方程进一步得到了广泛而有效地应用。如今，“物理问题变得越来越复杂，积分方程变得越来越有用”。 积分方程与数学的其他分支，例如，微分方程、泛函分析、复分析、计算数学、位势理论和随机分析等都有着紧密而重要地联系。甚至它的形成和发展是很多重要数学思想和概念的最初来源和模型。例如，对泛函分析中平方可积函数、平均收敛、算子等的形成，对一般线性算子理论的创立，以至於对整个泛函分析的形成都起着重要的推动作用。积分方程论中许多思想和方法，例如，关於第二种弗雷德霍姆(Fredholm)积分方程的弗雷德霍姆理论和奇异积分方程的诺特(Noether)理论以及逐次逼近方法，本身就是数学中经典而优美的理论和方法之一。 编辑本段起源 积分号下含有未知函数的方程。其中未知函数以线性形式出现的，称为线性积分方程；否则称为非线性积分方程。积分方程起源于物理问题。牛顿第二运动定律的出现，促进了微分方程理论的迅速发展，然而对积分方程理论发展的影响却非如此。1823年，N.H.阿贝尔在研究地球引力场中的一个质点下落轨迹问题时提出的一个方程，后人称之为阿贝尔方程，是历史上出现最早的积分方程，但是在较长的时期未引起人们的注意。“积分方程”一词是 P.du B.雷蒙德于1888年首先提出的。19世纪的最后两年，瑞典数学家（E.）I.弗雷德霍姆和意大利数学家V.沃尔泰拉开创了研究线性积分方程理论的先河。从此，积分方程理论逐渐发展成为数学的一个分支。 1899年，弗雷德霍姆在给他的老师（M.）G.米塔－列夫勒的信中，提出如下的方程 公式 , (1) 式中φ（x）是未知函数；λ是参数，K（x，y）是在区域0 ≤x，y≤1上连续的已知函数；ψ(x)是在区间0≤x≤1上连续的已知函数。并认为方程（1）的解可表为关于λ的两个整函数之商。1900年，弗雷德霍姆在 包括基本概念，差分格式的构造、截断误差和稳定性，这些内容是贯穿整个教材的主线。解答问题关键在过程，能够显示出你已经掌握了书上的内容，知道了解题方法。这次考试题目的类型：20分的选择题，主要是基本概念的理解，后面有五个大题，包括差分格式的构造、截断误差和稳定性。 习题一 1． 略 2． y y x f -=),(，梯形公式：n n n n n n y h h y y y h y y )121(),(2111+-+=+- =+++，所以0122)1(01])121[()121()121(y h h y h h y h h y h h n h h n n n +--+--+-+=+-+==+-+= ，当0→h 时， x n e y -→。 同理可以证明预报-校正法收敛到微分方程的解. 3． 局部截断误差的推导同欧拉公式; 整体截断误差: ? ++++++-++≤1 ),())(,(11111n n x x n n n n n n n dx y x f x y x f R εε 11)(++-++≤n n n y x y Lh R ε，这里R R n ≤ 而111)(+++-=n n n y x y ε，所以 R Lh n n += -+εε1)1(，不妨设1 6 数值积分 如果函数)(x f 在区间],[b a 上连续，且原函数为)(x F ，则可用牛 顿―莱布尼兹公式：)()()(a F b F dx x f b a -=?计算定积分。然而很多函数 无法用牛顿―莱布尼兹公式求定积分。 一个简单被积函数，例如错误!未找到引用源。dx cx bx a ?++2，其不定积分可能很复杂，见下面的MA TLAB 实例： >> syms a b c x >> int(sqrt(a+b*x+c*x*x),x) ans=1/4*(2*c*x+b)/c*(a+b*x+c*x^2)^(1/2)+1/2/c^(1/2)*log((1/2*b+c*x )/c^(1/2)+(a+b*x+c*x^2)^(1/2))*a-1/8/c^(3/2)*log((1/2*b+c*x)/c^(1/2)+(a+b*x+c*x^2)^(1/2))*b^2 所以有必要研究简单、高效的计算定积分的方法（即数值积分方法）。数值积分的基本思想是构造一个简单函数)(x P n 来近似代替被积分函数)(x f ，然后通过求?b a n dx x P )(得?b a dx x f )(的近似值。 6.1 插值型求积公式 设?=b a dx x f I )(* ，插值型求积公式就是构造插值多项式)(x P n ，使 ?=≈b a n dx x P I I )(*。 构造以a ，b 为结点的线性插值多项式)()()(1b f a b a x a f b a b x x P --+--= ，[])()()(21)()()(1b f a f a b dx b f a b a x a f b a b x dx x P T b a b a +-=?? ? ???--+--==??称为梯形公式。 微分方程数值解―― 第二章 习题 1． 设)('x f 为)(x f 的一阶广义导数，试用类似办法定义)(x f 的k 阶广义导数) () (x f k （ ,2,1=k ）。 解：对一维情形，函数的广义导数是通过分部积分来定义的。 我们知，)(x f 的一阶广义导数位)(x g ，如果满足 dx x x f dx x x g b a b a )()()()('?? -=?? 类似的，)(x f 的k 阶广义导数为)()() (x f x g k =，如果有 dx x x f dx x x g b a k k b a )()()1()()()(?? -=?? 2． 试建立与边值问题 ?????====<<=+=) 2.1(0)()(,0)()() 1.1(,''44b u b u a u a u b x a f u dx u d Lu 等价的变分问题。 证明： 设}0)()(,0)()(),(|{' '2====∈=b v a v b v a v I H v v V 对方程)1.1(两边同乘以v ，再关于x 在),(b a 上积分)(V v ∈，得 ??=+b a b a fvdx vdx u dx u d )(44 其中 dx dx dv dx u d dx dx dv dx u d dx u d v dx u d d v vdx dx u d b a b a b a b a b a ???? -=-==33 33333344|)( dx dx v d dx u d dx dv dx u d dx u d d dx dv b a b a b a ??+-=-=22222222|)( dx dx v d dx u d b a ? = 2 222 (*) 记dx uv dx v d dx u d v u a b a ?+=)(),(2 222，?=b a fvdx v f ),(。于是我们得到以下等价变分问题的提法： (1)：方程： 取x的范围为[0,10),a=0,b=5; 求解步骤： 1、对s进行离散，取,将[0,5]中每隔0.01取的数以此记为（i=0,1,2,3….500）,即有 2、对x进行离散，也取,将[0,10)中每隔0.01取的数以此记为（i=0,1,2,3….999）对于所取的任意有 写成矩阵的形式为 对于全部的，可写成： 1000 说明：s的范围可以为x的取值范围，也可以比x的取值范围小，当s的取值范围比x的取值范围小时，可以像上式那样将等式右边第二项的系数矩阵中s取不到的值令k=0，s和x的范围相等仅仅是一种特殊情况。 3、上式可化成F-KF=G的形式，进一步化成AF=G的形式，对其用doolittle分解法求解，可等到一系列离散值。 求解过程中取a=0,b=5,x的范围为[0,10)，， ，的理论值应为; C++程序如下： #include 实验09 数值微积分与方程数值求解 （第6章 MATLAB 数值计算） 一、实验目的 1. 掌握求数值导数和数值积分的方法。 2. 掌握代数方程数值求解的方法。 3. 掌握常微分方程数值求解的方法。 二、实验内容 1. 求函数在指定点的数值导数 232()1 23,1,2,302 6x x x f x x x x x == 程序及运行结果： 2. 用数值方法求定积分 (1) 22210 cos 4sin(2)1I t t dt π = ++? 的近似值。 程序及运行结果： 《数学软件》课内实验 王平 (2) 222 1I dx x π =+? 程序及运行结果： 3. 分别用3种不同的数值方法解线性方程组 6525494133422139211 x y z u x y z u x y z u x y u +-+=-??-+-=? ? ++-=??-+=? 程序及运行结果： 4. 求非齐次线性方程组的通解 123412341 2342736352249472 x x x x x x x x x x x x +++=?? +++=??+++=? 5. 求代数方程的数值解 (1) 3x +sin x -e x =0在x 0=1.5附近的根。 程序及运行结果（提示：要用教材中的函数程序line_solution ）： (2) 在给定的初值x 0=1，y 0=1，z 0=1下，求方程组的数值解。 23 sin ln 70 3210 50y x y z x z x y z ?++-=?+- +=??++-=? 6. 求函数在指定区间的极值 (1) 3cos log ()x x x x x f x e ++=在(0,1)内的最小值。 (2) 332 12112122(,)2410f x x x x x x x x =+-+在[0,0]附近的最小值点和最小值。 7. 求微分方程的数值解，并绘制解的曲线 22 50(0)0 '(0)0xd y dy y dx dx y y ?-+=??? =??=??? 程序及运行结果（注意：参数中不能取0，用足够小的正数代替）： 令y 2=y,y 1=y '，将二阶方程转化为一阶方程组： '112' 21 12 5 1(0)0,(0)0 y y y x x y y y y ?=-??=??==?? 8. 求微分方程组的数值解，并绘制解的曲线 123213 312123'''0.51(0)0,(0)1,(0)1 y y y y y y y y y y y y =??=-?? =-??===? 程序及运行结果：matlab数值微积分与方程数值求解
几种定积分的数值计算方法
(精选)实验二 数值方法计算积分
微分方程数值解法
偏微分方程数值解法答案
1. 积分方程一般概念与弗雷德霍姆方程
求第一类Fredholm积分方程的离散正则化方法
数值积分 (论文)
微分方程的分类及其数值解法
常微分方程数值解法
非线性Volterra积分方程(学习资料)
第一章积分方程的来源及基本概念
计算方法讲义：七 数值积分
积分方程
微分方程数值解法答案
工程中的计算方法课件6 数值积分
微分方程数值解――
积分方程的计算
实验09 数值微积分与方程数值解(第6章)
相关主题
文本预览