资源描述
计算方法（ A）大作业 姓名： 班级： 专业： 学号： 1 共轭梯度法 一、 算法原理 共轭梯度法是把求解线性方程组的问题转化为求解一个与之等价的二次函数极 小化的问题，因此从任意给定的初始点出发，沿一组关于矩阵 A 的共轭方向进行 线性搜索，在无舍入无差的假定下，最多迭代 n(其中 n 为矩阵 A 的阶数）次就可 求得二次函数的极小点，也就求得了线性方程组 Ax=B 的解。 下述定理给出了求系数矩阵 A 是对称正定矩阵的线性方程组 Ax=b 的解与求二 次函数 () = 12 极小点的等价性。 定理 3.4.1 设 A 是 n 阶对称正定矩阵，则 是方程组 Ax=b 的解的充分必要条 件是 是二次函数 () = 12 的极小点，即 = () = () 证明：充分性 .设 是 f(x)的极小点，则由无约束最优化问题最优解的必要条件知 () = 即 是方程组 Ax=b 的解。 必要性 . 若 是方程组 Ax=b 的解，即 A=b，注意到 A 是对称正定矩阵，故 x 有 ()() = 12 12 + = 12( 2 +) + = 12( 2() +)( ) = 12( )( ) 0 即 是 f(x)的极小点，进而由 A 是正定矩阵知， 是 f(x)的最小点。证毕。 共轭梯度法在形式上具有迭代法的特征，在给定初始值情况下，根据迭代公式： (+1) = () +() 产生的迭代序列 (1),(2),(3) 在无舍入误差假定下，最多经过 n 次迭代，就可 求得 f(x) 的最小值，也就是方程 Ax=b 的解。 共轭梯度法中关键的两点是，确定迭代格式中的搜索向量 ()和最佳步长 （ 0） 。 实际上，搜索方向 d()是关于矩阵 A 的共轭向量，在迭代中逐步构造之。 步长 的确定原则是给定迭代点 ()和搜索方向 d()后，要求选取非负实数 ，使 得 f() +d()达到最小，即选择 ，满足 () +() = 0 () +() 2 下面首先导出最佳步长 k 的计算公式。 假 设 迭 代 点 () 和 搜 索 方 向 d() 已 经 给 定 ， 可 以 通 过 一 元 函 数 () = () +() 的极小化 () = 0 () +() 来求得，根据多元复合函数的求导法则得 : () = () +()() 令 0 = () = () +()() = () +()() = () +()() = () +()() 其中 () = ()。由此解得最佳步长 = ()()()() (3.4.4) 下面确定搜索方向 d()。 给定初始向量 (0)后，由于负梯度方向是函数下降最快的方向，故第 1 次迭代 取搜索方向 (0) = (0) = (0) = (0) 。令 (1) = (0) +0(0) 其中 0 = (0)(0)(0)(0)。第 2次迭代时，从 (1) 出发的搜索方向不再取 (1)，而是 选取 (1) = (1) +0(0)，使得 (1)与 (0)是关于矩阵 A 的共轭向量，即要求 (1)满 足 (1),(0) = 0，由此可求得参数 0: 0 = (1)(0) (0)(0) 然后从 (1)出发，沿 (1)进行搜索得到 (2) = (1) +1(1) 其中 1由式 （ 3.4.4） 确定。 一般地， 设已经求出 (+1) = () +()，计算 (+1) = (+1)。令 (+1) = (+1) +()，选取 ，使得 (+1)和 ()是关于 A 的共轭向量，即要求 (+1)满足 (+1),() = 0。由此可得 : = (+1)() ()() 从而确定了搜索方向 (+1)。 3 综上所述，共轭梯度法的计算公式为 (0) = (0) = (0), = ()()()(), (+1) = () +(), (+1) = (+1), = (+1)()()() , (+1) = (+1) +(). 利用定理 3.4.3，可以将 和 的表达式简化为 = ()() ()() = (+1)(+1) ()() 定理 3.4.3 设分别是共轭梯度法中产生的非零残向量和搜索方向，则 （ 1） (),(1) = 0 （ 2） (),() = (),() （ 3） (0), (1), ()(k n1)是正交向量组， (0), (1), ()(k n1) 是关于 A的共轭向量组。 关于共轭向量法的误差有如下定理： 定理 3.4.5 设是用共轭梯度法求得的迭代序列，则有误差估计式 () 2( 1 +1) (0) 其中范数 = , = 2(). 若 () = 0，则 ()就是线性方程组的解，故在计算过程中将 (+1) 作为 迭代终止的条件。 共轭梯度法的 计算过程如下： （ 1） 给定初始近似向量 (0) 以及精度 0； （ 2） 计算 (0) = (0)，取 (0) = (0)； （ 3） for k=0 to n-1 do （ i） = ()()()()； （ ii） (+1) = () +()； （ iii） (+1) = (+1)； （ iv） 若 (+1) 或 k+1=n，则输出近似解 (+1)，停止；否则，转（ v）； 4 （ v） = (+1) 2 2 ()22 ； （ vi） (+1) = (+1) +()； end do 计算经验表明，对于不是十分病态的问题，共轭梯度法收敛较快，迭代次数远 小于矩阵的阶数 n。对于病态问题，只要进行足够多次的迭代（迭代次数大约为矩 阵阶数 n 的 3-5 倍）后，一般也能得到满意的结果，因而共轭梯度法是求解高阶 稀疏线性方程组的一个有效、常用的方法 。 二、 程序框图 5 三、程序使用说明及案例计算结果 程序使用说明 本共轭梯度法求解线性方程的程序需要用户给定对称正定矩阵 A 的阶数，误 差以及初始向量，然后程序自动调用定义的函数 Gongetidu(A,b,x0,epsa)实现求解线 性方程组 Ax=b。其中 A， b 的阶数，初始向量 x0 和 epsa 都是可变的。 案例计算结果 可以发现，当 n=100， 200， 300 时，方程组 Ax=b 的解 x 为元素均为 1 的 n 阶向量。 6 四、 Matlab 程序（ M-文件） clear A b x0;clc; %程序运行前首先清除 A， b 和 X0 的数值，以免影响计算 n=input(请输入对称正定矩阵 A 的阶数 n=); epsa=input(请输入误差 =); A=zeros(0,0); %给矩阵 A 预先配置内存空间，减少耗时 for k=1:(n-1) %读取题目 3.2 的矩阵 A A(k,k)=-2; A(k,k+1)=1; A(k+1,k)=1; end A(n,n)=-2; A; b(1,1)=-1; %读取题目 3.2 的矩阵 b b(n,1)=-1; b; x0(1,1)=input(假定初始向量各元素相同，均为： ); %给题目 3.2 的迭代过程赋初值 for kk=2:n x0(kk,1)=x0(1,1); end x=Gongetidu(A,b,x0,epsa); %调用共轭梯度法求线性方程组的函数 function x=Gongetidu(A,b,x0,epsa) n=size(A,1); x=x0; r=b-A*x; d=r; for k=0:(n-1) alpha=(r*r)/(d*A*d); x=x+alpha*d; r2=b-A*x; if (norm(r2)=epsa)|(k=n-1) disp(方程组的近似解为 x=); disp(x); break; end beta=norm(r2)2/norm(r)2; d=r2+beta*d; r=r2; end end 7 三次样条插值 三、 算法原理 分段低次插值多项式的一阶或二阶导数不连续，不能满足许多实际问题的要求。 例如，在飞机机翼外形设计，船体放样等许多实际问题中，要求曲线二阶光滑，即 要求函数的二阶导数连续。为了将一些已知点连成一条充分光滑的曲线，早期的做 法是绘图人员用压铁把一条称为样条的富有弹性的细木条或金属条固定在相邻的若 干点上，然后沿样条画下一条所需的光滑曲线，这条曲线称为样条曲线。实际上， 它是由分段三次曲线拼接而成的，在连接点处二阶导数连续。将样条曲线抽象为数 学问题称为样条函数。 1、三次样条插值函数的定义 设在区间 a,b上给定 n+1 个 节点 ( 0 1 )，在节点 处的 函数值 = ()( = 0,1,)。若函数 S(x)满足以下三条： （ 1） 在每个子区间 1,上， S(x)是三次多项式； （ 2） 是个 S() = ( = 0,1,)； （ 3） 在区间 a,b上， S(x)的二阶导数 () 连续， 则称 S(x)为函数 y = f(x)在区间 a,b上的三次样条插值函数。 由定义可知， S(x)是分段三次多项式，故在每个子区间 1, ( = 1,)上， S(x)有 4 个待定参数，共有 n 个子区间，所以 S(x)共有 4n 个待定参数。根据定义中 的条件（ 3）知 () = +()， () = + ()， = 1,2, 1. () = + ()， 上式共有 3n-3 个条件，再加上定义条件（ 2）中的 n+1 个条件，共有 4n-2 个 条件。还需要增加两个条件才能确定 4n 个待定参数，即才能确定 S(x)。所增加的 条件称为边界条件或端点条件。三种常用的边界条件如下： （ 1） 已知 f(x)在区间 a,b的两端点 a， b 处的二阶导数值 () ,() ； （ 2） 已知 f(x)在区间 a,b的两端点 a， b 处的一阶导数值 () ,() ； （ 3） 已知函数 f(x)式以 T=b-a 为周期的周期函数。 2、三次样条插值函数的导出 1) 导出在子区间 ,( = ,) 上的 S(x)的表达式 由于 S(x)的二阶导数连续，设 S(x)在 节点 处的二阶导数值 () = ( = 1,2,)，其中 为未知的待定参数。由 S(x)是分段的三次多项式知， () 是分段 线性函数， () 在子区间 1,上可表示为 8 () = 1 1 + 1 1 = 1 + 1 1 其中 = 1。 对上式 两端 积分 两次 得 S(x) = ( ) 3 6 1 + ( 1)3 6 +( )+( 1) 1 其中 ,为积分常数。现在确定 ,。 由插值条件 S(1) = 1,S() = 得 2 6 1 + = 1, 2 6 + = . 由此解得 = (1 2 6 1)/ = ( 2 6 )/ 将 , 代入 S(x)可得 S(x)在子区间 1,( = 1,2,)上的表达式为 S(x) = ( ) 3 6 1 + ( 1)3 6 +(1 2 6 1) +( 2 6 ) 1 1 2)建立参数 的方程组 对 S(x)求导可得 S(x) = ( ) 2 2 1 + ( 1)2 2 + 1 16 1 上式中令 = 得 S(x)在 处的左导数 S () = 2 + 1 16 = 6 1 +3 + 1 令 = 1得 S(x)在 1处 的 右导数 S+ (1) = 2 1 + 1 16 = 3 1 6 + 1 9 从而 S+ () = +13 +16 +1 +1 +1 由 S(x)在内节点 处一阶导数 的 连续 性可知 () = + () = 1,2, 1. 即 6 1 + +1 3 + +1 6 +1 = +1 +1 1 = 1,2, 1. 上式两端同乘 6 +1 得 +1 1 +2 + +1 +1 +1 = 6 +1 ( +1 +1 1 ) 记 = +1 = +1 +1 = 1 = 6 +1 (+1 +1 1 ) = 61,+1 = 1,2, 1. 则关于 的方程组可写为 1 +2 +1 = = 1,2, 1. 上式 中的未知量 在力学上解释为细梁在 截面处的弯矩，故称为三弯矩方程 组。 三弯矩方程组 中 只有 n-1 个方程， 不能 完全 确定 n+1 个未知量 ( = 0,1,2,)，所以 欲解三弯矩方程组或者减少两个未知量，或者再增加两个方程， 这由 边界条件 来 确定。 3)三种边界条件的三弯矩方程 (1)第一种边界条件 ： 已知 f(a) ,f(b) 。 取 0 = f(a) , = f(b) ，这时三弯矩方程组 减少了两个未知量，变成只含 n-1 个未知量 ( = 1,2, 1)的 n-1 个方程的方 程组 21 +12 = 1 00 1 +2 +1 = = 2,3, 2. (1) 12 +21 = 1 1 用矩阵表示为 ( 2 1 2 2 2 2 2 2 1 2 ) ( 1 2 2 1) = ( 1 10 2 2 1 1) 10 注意到 + = 1 2，以上方程组系数矩阵是 严格三对角占优矩阵，可用追 赶法求解。 解出 ( = 1,2, 1)后，代入 S(x)即 可得三次样条插值函数的数学 表达式。 (2)第二种边界条件：已知 f(a) ,f(b) 。记 0 = f(a) , = f(b) ，则有 S+ (0) = 0 ,S () = 。故有 13 0 16 1 +1 0 1 = 0 6 1 +3 + 1 = 即 20 +1 = 0 1 +2 = 其中 0 = 6 1 (1 0 1 0 ) = 60,0,1 = 6 ( 1 ) = 61, 将以上两个方程添加稻方程组中，得到关于第二种边界条件的三弯矩方程组 20 +1 = 0 1 +2 +1 = = 1,2,3, 1. (2) 1 +2 = 用矩阵表示为 ( 2 1 1 2 1 2 2 2 1 2 1 1 2 ) ( 0 1 2 1 ) = ( 0 1 2 1 ) 同样，该方程的系数矩阵是严格三对角占优矩阵，可用追赶法求解。 (3)第三种边界条件：周期型边界条件 。 已知 y=f(x)是以 T = ba = 0 为周期 的周期函数，则由周期性可知 = 0,+1 = 1, = 0,+1 = 1,+1 = 1。 这时， 将点 看成内 节 点，则方程组对 i=n 也成立，即 有 1 +2 +1 = 也即 1 +1 +2 = 其中 = 6 +1 (+1 +1 1 ) = 6 +1 (1 1 1 ) 方程组的第 1 个（ i=1）方程为 10 +21 +12 = 1 11 即 21 +12 +1 = 1 于是 可得关于第三种边界条件的三弯矩方程组 21 +12 +1 = 1 1 +2 +1 = = 2,3, 1, (3) 1 +1 +2 = 用矩阵表示为 ( 2 1 1 2 2 2 1 2 1 2 ) ( 1 2 1 ) = ( 1 2 1 ) 显然，该方程的系数矩阵也是严格对角占优矩阵，可用 LU 方法求解。 根据不同的边界条件，求解相应的三弯矩方程组，将求得的 代入 S(x)的表达 式便可得三次样条插值函数 S(x)。 四、 程序框图 12 三、程序使用说明及案例计算结果 改程序不需要用户输入数据，可直接计算 5()、 10()、 10()，并将其分 别与原函数 () = 11+252画在同一坐标系内，可以清楚的观察到 5()的拟合效 果没有 10()好，而 10()的和原函数基本重合，误差最小。 程序运行结果如下： 13 四、 Matlab 程序（ M-文件） 改程序分为主程序、 threetimes 函数、 newtonway 函数三部分，分别如下： 主程序 clear;clc; threetimes(10) %调用三次样条插值 段数为 10 hold on newtonway(5) %调用牛顿插值法 段数为 5 hold on newtonway(10) %调用牛顿插值法 段数为 10 hold on threetimes 函数 function threetimes(n) k=1;syms t; y=1/(1+25*k.2); %测试函数公式 %区间段数，由函数传递值确定 h=2/n; %区间段长 cout=1; %cout 为计数器 x=; %x 轴节点数 for i=-1:h:1 %计算节点数值 k=i; x(cout)=i; y(cout)=1/(1+25*k2); 14 cout=cout+1; end % x=-3 -1 0 3 4; % y=7 11 26 56 29; for i=1:1:size(x,2) %b 插商矩阵 前两列赋值 b(i,1)=x(i); b(i,2)=y(i); end f=1; %计数变量 for j=3:1:size(x,2)+1 %计算插商循环 size(x,2)区间节点个数 k=1; for i=j-1:1:size(x,2) b(i,j)=(b(i,j-1)-b(i-1,j-1)/(b(i,1)-b(i+2-j,1); if(k=1) c(f)=b(i,j);%提取插商系数 k=0; f=f+1; end end end for i=3:1:size(x,2) %d 等于 6 倍于 三阶插商 d(i-2)=6*b(i,4); %d 为求解 m 矩阵的最后一侧列向量 end for i=2:1:size(x,2) h(i-1)=b(i,1)-b(i-1,1); end for i=1:1:size(h,2)-1 u(i)=h(i)./(h(i)+h(i+1); r(i)=1-u(i); end for i=1:1:size(x,2)-3 bb(i,i)=2; %对角阵数值 bb(i,i+1)=r(i); %对角数上面的值 bb(i+1,i)=u(i+1); %对角下面的数值 end bb(size(x,2)-2,size(x,2)-2)=2; mm=bbd; m=; % M 矩阵 M（ 0）和 M（ n） =0; m(1)=0; 15 for i=1:1:size(mm,1) m(i+1)=mm(i); end m(size(m,2)+1)=0; subplot() for i=2:1:size(x,2) % 通过公式确定每段拟合函数 s_1=m(i-1)*(x(i)-t).3/(6*h(i-1); % 注意其中下标顺序 s_2=m(i)*(t-x(i-1).3/(6*h(i-1); s_3=(y(i-1)-h(i-1).2*m(i-1)/6)*(x(i)-t)/h(i-1); s_4=(y(i)-h(i-1).2*m(i)/6)*(t-x(i-1)/h(i-1); s=s_1+s_2+s_3+s_4; s=expand(s); a=x(i-1):0.01:x(i); f=subs(s,t,a); plot(a,f,r-.,linewidth,3) hold on end x=-1:0.01:1; axis(-1 1 0 1.1); y1=power(1+25*power(x,2),-1); plot(x,y1,b-,linewidth,1) xlabel(Variableitx) ylabel(Variableity) grid hold on end newtonway 函数 function newtonway(n) k=1;syms t; y=1/(1+25*k.2); %公式 %区间段数 h=2/n; %区间短长 cout=1; x=; for i=-1:h:1 k=i; x(cout)=i; y(cout)=1/(1+25*k2); %计算节点数值 cout=cout+1; 16 end size(x,2); b=size(x,2),n; for i=1:1:size(x,2) %b 插商矩阵赋值 b(i,1)=x(i); b(i,2)=y(i); end f=1; %计数变量 for j=3:1:size(x,2)+1 %计算插商循环 m=1; for i=j-1:1:n+1 b(i,j)=(b(i,j-1)-b(i-1,j-1)/(b(i,1)-b(i+2-j,1); if(m=1) c(f)=b(i,j);%提取插商系数 m=0; f=f+1; end end end b; l=1;d=b(1,2);d; for i=1:1:size(x,2)-1 %计算出多项式 l=l*(t-b(i,1); d=d+c(i)*l; end d=simplify(d); s=2; %设定画图区间上下界数值 t=-s:0.01:s; k=-s:0.01:s; y=1./(1+25.*k.*k); figure; axis(-1,1,-0.5,2); set(gca,XTickMode,manual,XTick,0 1 2 3 4 5 6 7 8 9 10); set(gca,YTickMode,manual,YTick,-1 0 1);grid xlabel(Variableitx) ylabel(Variableity) hold on plot(t,subs(d),r-.,k,y,b) end 17 龙格贝积分 五、 算法原理 对于复化梯形求积公式而言，可取积分的近似值为 I(f) 2 + 141(2 ) = 2 (1) 对于复化辛普森求积公式有 I(f) = 2880 4(4)() I(f)2 = 2880 (2)4(4)(1) 1 假定 (4)()在区间 a,b上连续且变化不大，可以认为 (4)() (4)(1)。将以上两 式相除得 I(f) I(f)2 16 由此解得 I(f) 2 + 142 1(2 ) 2 (2) 同理，对复化科茨求积公式有 I(f) 2 + 143 1(2 ) 2 (3) 下面对 (1)、 (2)、 (3)作进一步分析，看一看它们的实质和它们之间的关系。 2 = 2 +13(2 ) = 13(42 ) = 13412 +2(1 +2 ) =1 = 13 +2(1 +2 ) =1 = 132()+2()+() 1 =1 +2(1 +2 ) =1 = 16()+2()+4(1 +2 ) =1 +() 1 =1 = 即 = 2 + 141(2 ) 同理可得 18 = 2 + 142 1(2 ) 令 = 2 + 143 1(2 ) (4) (4)式即称为龙格贝积分公式。其截断误差为 8(8)()。以此类推可知，若取 2 + 144 1(2 ) 作为积分的近似值其截断误差是 c10(10)()。容易设想取其作为积分的近似值可 能更精确。事实上并非完全如此，这是因为被积函数 f(x)的高阶导数性态很难估计， 特别是高阶导数性态不好的函数， (10)()的数量级可能很大，虽然 10很小，但截 断误差项 c10(10)()可能很大。如此计算的效果可能并不理想。再者，一般地， 2 是两个相近的数相减，可能产生较大的误差。因此，一般不再继续进行外 推，计算到龙格贝积分公式就已经足够准确。如果积分区间 a,b较大，则最好事先 把 a,b分成若干个子区间，然后在每个子区间上再进行数值积分。 龙格贝积分法是将区间 a,b逐次分半，逐次递推计算，再用 (1)(4)逐次计算， 以得到较高精度的积分近似值。龙格贝积分法的计算公式如下： 1 = 2 ()+() 2+1 = 122 + 2+1 ( +(2 1) 2+1 ) 2 =1 = 0,1,2, 2 = 2+1 + 141(2+1 2) 2 = 2+1 + 142 1(2+1 2) 2 = 2+1 + 143 1(2+1 2) 计算过程如下图所示进行。 1 2 1 22 2 1 23 22 2 1 24 23 22 2 19 若 |2+1 2| ，则取 If 2+1；否则继续计算， 直到 满足精度要求为止。 六、 程序框图 开 始 输 入 , , , b 3k 1kk fa 23 22 kkRR T 输 出 近 似 解 2 2 kI f R F 结 束 1 ( ) ( ) 2 ba T f a f b 1 2 1122 1 1 ( ( 2 1 ) ) , 0 , 1 , 2 2 2 2 k kk kk i b a b a T T f a i k 11 2 2 2 2 1 ( ) , 0 , 1 , 2 41 k k k kS T T T k 11 22 2 2 2 1 ( ) , 0 , 1 41 k k k kC S S S k 11 32 2 2 2 1 ( ) , 0 41 k k k kR C C C k 1 2 1122 1 1 ( ( 2 1 ) ) 2 2 2 k kk kk i b a b a T T f a i 11 2 2 2 2 1 () 41 k k k kS T T T 11 22 2 2 2 1 () 41 k k k kC S S S 2 1 1 2 32 2 2 2 1 () 41 k k k kR C C C 三、程序使用说明及案例计算结果 改程序只需要输入精度要求即可求解习题 6.2 中的 4 个小题，其中第 3、第 4 个小题的函数边界值比较特殊，故进行近似化处理。程序运行结果如下： 20 四、 Matlab 程序（ M-文件） clear;clc; esp=input(请输入精度 esp=); for i=1:4 if i=eps R(2k)=C(2(k+1),a,b,n)+1/(43-1)*(C(2(k+1),a,b,n)-C(2k,a,b,n); eps1=abs(R(2k)-R(2(k-1); k=k+1; end y=R(2(k-1); end 22 四阶龙格 -库塔法 七、 算法原理 利用泰勒展开可以导出龙格 -库塔法（简称 R-K 方法）。 m 级 的龙格 -库塔法的 一般形式为 +1 = +11 +22 + 1 = (,) 2 = ( +2, +211) 3 = ( +3, +311 +322) = ( +, +11 +22 +,11) 其中 ,均为常数，有待定系数法确定。确定的原则是将局部截断误差 = (+1)+1在处 泰勒展开，适当选取 h 的系数，使得局部截断误差 的阶数 尽可能高。 下面以 m=2 级的 R-K 法为例来说明 R-K 法的导出思想。 2 级 R-K 法形式如下： +1 = +11 +22 1 = (,) 2 = ( +, +1) 将 (+1)在 处泰勒展开有 (+1) = ( +) = ()+() + 12!()2 + 13!()3 +(4) 由 () = (,()知 () = + = + () = +2 +2 + + 2 所以 (+1) = ()+ +12( + )2 +16( +2 +2 + + )3 +(4) 其中 为 (,()的简写，符号 , ,的含义亦然。 又由二元函数的泰勒展开式有 2 = ( +, +1) = ( + + +1222 +2 +12222 +) 于是 +1 = +11 +22 23 = +1 +2( + + +1222 +2 +12222 +) = +(1 +2) +2( + )2 +2 (122 + +1222)3 +(4) 可得局部截断误差 = (+1)+1 = (11 2) +(122) +(122) 2 +(161222) +(132) +(161222)2 +16( + 2)3 +(4) 为使 的阶尽可能高，由 的任意性，应选取 ,1,2，使得 ,2的系数为零， 即得方程组 11 2 = 0 1 22 = 0 1 22 = 0 由此解出 ,1,2，代入相应表达式即可。注意到 ,1,2无论怎样选取， 3的 系数中 + 2 0，故 = (3)，所以通过解上面的方程组得到的方法均 是二阶方法。其含有 3 个方程， 4 个未知量，解有无穷个。 (1)取 1 = 2 = 12, = = 1，则得 +1 = +121 +122 1 = (,) 2 = ( +, +1) 即改进欧拉法。 (2)取 1 = 0,2 = 1, = = 1/2，则可得 +1 = +2 1 = (,) 2 = ( +12, +121) 通常也写为 +1 = +( +12, +12(,) 即变形欧拉法。 24 类似于上述二阶方法的推导，可得多种 4 级 4 阶 R-K 法。应用最广泛的是如下 标准（经典） 4 级 4 阶 R-K 法： +1 = +16(1 +22 +23 +4) 1 = (,) 2 = ( +12, +121) 3 = ( +12, +122) 4 = ( +, +3) 另外两个常用的 4 级 4 阶 R-K 法的计算公式不再列出。 4 阶以上的 R-K 法计算函数值的工作量较大，对于大量的实际问题 4 阶 R-K 法 已可满足精度要求。若有更高的精度要求，可使用 6 级 5 阶英格兰公式（此处不 再给出）。 R-K 方法的优缺点如下： 优点：精度高，它是显式单步法，计算 +1只需要 一个值，每一步的步长可根据 精度的要求单独考虑其大小。 缺点：需要计算多个函数值，当 (,)比较复杂时，计算量较大。 八、 程序框图 原题目是高阶微分方程初值问题，需将其转化为一阶微分方程组来求解。令 1 = ,2 = ,3 = ,转化后的方程如下： 1 = 2 2 = 3 3 = 3 +2 1 +2 3 1(0) = 1,2(0) = 3,3(0) = 3 1,2,3,4计算如下： 1 = (,) = ( 11 12 13 ) 2 = ( +12, +121) = ( 21 22 23 ) 3 = ( +12, +122) = ( 31 32 33 ) 4 = ( +, +3) = ( 41 42 43 ) 25 程序框图如下： 开始 输入区间 a,b，步长 h 初值 1(),2(),3() 1 = (,) = ( 11 12 13 ) 2 = ( +12, +121) = ( 21 22 23 ) 3 = ( +12,
展开阅读全文