第8章 MATLAB方程求解

第8章MATLAB方程求解

8.1线性方程组求解8.2非线性方程数值求解8.3最优化问题求解8.4常微分方程初值问题的数值求解第8章MATLAB方程求解8.1线性方程组求解在MATLAB中，关于线性方程组的解法一般分为两类：一类是直接法，就是在没有舍入误差的情况下，通过有限步的矩阵初等运算来求得方程组的解；另一类是迭代法，就是先给定一个解的初始值，然后按照一定的迭代算法进行逐步逼近，求出更精确的近似解。第8章MATLAB方程求解8.1.1线性方程组的直接解法在MATLAB中，这些算法已经被编成了现成的库函数或运算符，因此，只需调用相应的函数或运算符即可完成线性方程组的求解。1．利用左除运算符的直接解法线性方程组求解最简单的方法就是使用左除运算符“\”，系统会自动根据输入的系数矩阵判断选用哪种方法进行求解。对于线性方程组Ax=b，可以利用左除运算符“\”求解：x=A\b。第8章MATLAB方程求解例8-1用直接解法求解下列线性方程组。A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];b=[13,-9,6,0]';x=A\b程序运行结果为：x=-66.555625.6667-18.7778

26.5556第8章MATLAB方程求解2．利用矩阵的分解求解线性方程组（1）LU分解MATLAB提供的lu函数用于对矩阵进行LU分解，其调用格式为：①[L,U]=lu(X)：产生一个上三角阵U和一个变换形式的下三角阵L（行交换），使之满足X=LU。注意，这里的矩阵X必须是方阵。②[L,U,P]=lu(X)：产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P，使之满足PX=LU。当然矩阵X同样必须是方阵。当使用第一种格式时，矩阵L往往不是一个下三角阵，但可以通过行交换成为一个下三角阵。第8章MATLAB方程求解设：则对矩阵A进行LU分解的命令如下：>>A=[1,-1,1;5,-4,3;2,1,1];>>[L,U]=lu(A)L=0.2000-0.07691.00001.0000000.40001.00000U=5.0000-4.00003.000002.6000-0.2000000.3846为检验结果是否正确，输入命令：>>LU=L*ULU=

1-115-43211第8章MATLAB方程求解利用第二种格式对矩阵A进行LU分解。>>[L,U,P]=lu(A)L=1.0000000.40001.000000.2000-0.07691.0000U=5.0000-4.00003.000002.6000-0.2000000.3846P=010001100>>LU=L*U %这种分解其乘积不为ALU=5-432111-11>>inv(P)*L*U %考虑矩阵P后的结果ans=1-115-43211实现LU分解后，线性方程组Ax=b的解x=U\(L\b)或x=U\(L\Pb)，这样可以大大提高运算速度。第8章MATLAB方程求解例8-2用LU分解求解例8-1中的线性方程组。A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];b=[13,-9,6,0]';[L,U]=lu(A);x=U\(L\b)程序运行结果为：x=-66.555625.6667-18.777826.5556或采用LU分解的第二种格式，命令如下：>>[L,U,P]=lu(A);>>x=U\(L\P*b)第8章MATLAB方程求解（2）QR分解对矩阵X进行QR分解，就是把X分解为一个正交矩阵Q和一个上三角阵R的乘积形式。QR分解只能对方阵进行。MATLAB的函数qr可用于对矩阵进行QR分解，其调用格式为：①[Q,R]=qr(X)：产生一个正交矩阵Q和一个上三角阵R，使之满足X=QR。②[Q,R,E]=qr(X)：产生一个正交矩阵Q、一个上三角阵R以及一个置换矩阵E，使之满足XE=QR。第8章MATLAB方程求解设：则对矩阵A进行QR分解的命令如下：>>A=[1,-1,1;5,-4,3;2,7,10];>>[Q,R]=qr(A)Q=-0.1826-0.0956-0.9785-0.9129-0.35320.2048-0.36510.9307-0.0228R=-5.47721.2780-6.572708.02298.151700-0.5917>>QR=Q*RQR=

1.0000-1.00001.00005.0000-4.00003.00002.00007.000010.0000第8章MATLAB方程求解利用第二种格式对矩阵A进行QR分解：>>[Q,R,E]=qr(A)Q=-0.0953-0.2514-0.9632-0.2860-0.91990.2684-0.95350.30110.0158R=-10.4881-5.4347-3.432506.0385-4.2485000.4105E=001010100>>Q*R/E%验证A=Q*R*inv(E)ans=1.0000-1.00001.00005.0000-4.00003.00002.00007.000010.0000实现QR分解后，线性方程组Ax=b的解x=R\(Q\b)或x=E(R\(Q\b))。第8章MATLAB方程求解例8-3用QR分解求解例8-1中的线性方程组。程序如下：A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];b=[13,-9,6,0]';[Q,R]=qr(A);x=R\(Q\b)程序运行结果为：x=-66.555625.6667-18.777826.5556或采用QR分解的第二种格式，命令如下：>>[Q,R,E]=qr(A);>>x=E*(R\(Q\b))将得到与上面同样的结果。第8章MATLAB方程求解（3）Cholesky分解如果矩阵X是对称正定的，则Cholesky分解将矩阵X分解成一个下三角阵和一个上三角阵的乘积。设上三角阵为R，则下三角阵为其转置，即X=R'R。MATLAB函数chol(X)用于对矩阵X进行Cholesky分解，其调用格式为：①R=chol(X)：产生一个上三角阵R，使R'R=X。若X为非对称正定，则输出一个出错信息。②[R,p]=chol(X)：这个命令格式将不输出出错信息。若X为对称正定的，则p=0，R与上述格式得到的结果相同，否则p为一个正整数。如果X为满秩矩阵，则R为一个阶数为q=p-1的上三角阵，且满足R‘R=X(1:q,1:q)。第8章MATLAB方程求解设：则对矩阵A进行Cholesky分解的命令如下：>>A=[2,1,1;1,2,-1;1,-1,3];>>R=chol(A)R=1.41420.70710.707101.2247-1.2247001.0000可以验证R'R=A，命令如下：>>R'*Rans=2.00001.00001.00001.00002.0000-1.00001.0000-1.00003.0000第8章MATLAB方程求解利用第二种格式对矩阵A进行Cholesky分解：>>[R,p]=chol(A)R=1.41420.70710.707101.2247-1.2247001.0000p=0结果中p=0，这表示矩阵A是一个正定矩阵。如果试图对一个非正定矩阵进行Cholesky分解，则将得出错误信息，所以，chol函数还可以用来判定矩阵是否为正定矩阵。实现Cholesky分解后，线性方程组Ax=b变成R'Rx=b，所以x=R\(R'\b)。第8章MATLAB方程求解例8-4用Cholesky分解求解例8-1中的线性方程组。命令如下：>>A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];>>b=[13,-9,6,0]';>>R=chol(A)错误使用chol矩阵必须为正定矩阵。命令执行时，出现错误信息，说明A为非正定矩阵。第8章MATLAB方程求解8.1.2线性方程组的迭代解法迭代法是一种不断用变量的旧值递推新值的过程，是用计算机解决问题的一种基本方法。它利用计算机运算速度快、适合做重复性操作的特点，让一组指令重复执行，在每次执行这组指令时，都从变量的原值推出它的新值。迭代解法非常适合求解大型稀疏矩阵的方程组。在数值分析中，迭代解法主要包括Jacobi迭代法、Gauss-Serdel迭代法、超松弛迭代法和两步迭代法。第8章MATLAB方程求解迭代法的思想:将方程改写为：这种形式的好处是将一组x代入右端，可以立即得到另一组x。如果两组x相等，那么它就是方程组的解，不等时可以继续迭代。可以构造方程的迭代公式为：第8章MATLAB方程求解1．Jacobi迭代法对于线性方程组Ax=b，如果A为非奇异方阵，且aii≠0（i=1，2，…，n），则可将A分解为A=D-L-U，其中D为对角阵，其元素为A的对角元素，L与U为A的下三角阵取反和上三角阵取反：于是Ax=b转化为：x=D-1(L+U)x+D-1b与之对应的迭代公式为：x(k+1)=D-1(L+U)x(k)+D-1b这就是Jacobi迭代公式。如果序列{x(k+1)}收敛于x，则x必是方程Ax=b的解。第8章MATLAB方程求解Jacobi迭代法的MATLAB函数文件jacobi.m如下：function[y,n]=jacobi(A,b,x0,ep)ifnargin==3ep=1.0e-6;elseifnargin<3errorreturnendD=diag(diag(A)); %求A的对角矩阵L=-tril(A,-1); %求A的下三角阵U=-triu(A,1); %求A的上三角阵B=D\(L+U);f=D\b;y=B*x0+f;n=1; %迭代次数whilenorm(y-x0)>=epx0=y;y=B*x0+f;n=n+1;end第8章MATLAB方程求解例8-5用Jacobi迭代法求解下列线性方程组。设迭代初值为0，迭代精度为10-6。在程序中调用函数文件jacobi.m，程序如下：A=[10,-1,0;-1,10,-2;0,-2,10];b=[9,7,6]';[x,n]=jacobi(A,b,[0,0,0]',1.0e-6)程序运行结果为：x=0.99580.95790.7916n=11第8章MATLAB方程求解2．Gauss-Serdel迭代法在Jacobi迭代过程中，计算

时，~已经得到，不必再用~

，即原来的迭代公式Dx(k+1)=(L+U)x(k)+b可以改进为Dx(k+1)=Lx(k+1)+Ux(k)+b于是得到：x(k+1)=(D-L)-1Ux(k)+(D-L)-1b该式即为Gauss-Serdel迭代公式。和Jacobi迭代相比，Gauss-Serdel迭代用新分量代替旧分量，精度会高些。第8章MATLAB方程求解Gauss-Serdel迭代法的MATLAB函数文件gauseidel.m如下：function[y,n]=gauseidel(A,b,x0,ep)ifnargin==3ep=1.0e-6;elseifnargin<3errorreturnendD=diag(diag(A)); %求A的对角矩阵L=-tril(A,-1); %求A的下三角阵U=-triu(A,1); %求A的上三角阵G=(D-L)\U;f=(D-L)\b;y=G*x0+f;n=1;%迭代次数whilenorm(y-x0)>=epx0=y;y=G*x0+f;n=n+1;end第8章MATLAB方程求解例8-6用Gauss-Serdel迭代法求解例8-5中的线性方程组。在程序中调用函数文件gauseidel.m，程序如下：A=[10,-1,0;-1,10,-2;0,-2,10];b=[9,7,6]';[x,n]=gauseidel(A,b,[0,0,0]',1.0e-6)程序运行结果为：x=0.99580.95790.7916n=7由此可见，一般情况下Gauss-Serdel迭代比Jacobi迭代要收敛快一些。但这也不是绝对的，在某些情况下，Jacobi迭代收敛而Gauss-Serdel迭代却可能不收敛。

第8章MATLAB方程求解例8-7分别用Jacobi迭代和Gauss-Serdel迭代法求解下列线性方程组，看是否收敛。>>a=[1,2,-2;1,1,1;2,2,1];>>b=[9;7;6];>>[x,n]=jacobi(a,b,[0;0;0])x=-27268n=4>>[x,n]=gauseidel(a,b,[0;0;0])x=NaNNaNNaNn=

1012可见对此方程，用Jacobi迭代收敛，而Gauss-Serdel迭代不收敛。因此，在使用迭代法时，要考虑算法的收敛性。第8章MATLAB方程求解8.1.3求线性方程组的通解①当系数矩阵A是一个满秩方阵时，方程Ax=b称为恰定方程，方程有唯一解x=A-1b，这是最基本的一种情况。一般用x=A\b求解速度更快。②当方程组右端向量b=0时，方程称为齐次方程组。齐次方程组总有零解，因此称解x=0为平凡解。当系数矩阵A的秩小于n（n为方程组中未知变量的个数）时，齐次方程组有无穷多个非平凡解，其通解中包含n-rank(A)个线性无关的解向量，用MATLAB的函数null(A,'r')可求得基础解系。③当方程组右端向量b≠0时，系数矩阵的秩rank(A)与其增广矩阵的秩rank([A,b])是判断其是否有解的基本条件：（a）当rank(A)=rank([A,b])=n时，方程组有唯一解：x=A\b或x=inv(A)*b。（b）当rank(A)=rank([A,b])<n时，方程组有无穷多个解，其通解=方程组的一个特解+对应的齐次方程组Ax=0的通解。可以用A\b求得方程组的一个特解，用null(A,'r')求得该方程组所对应的齐次方程组的基础解系，基础解系中包含n-rank(A)个线性无关的解向量。（c）当rank(A)<rank([A,b])时，方程组无解。第8章MATLAB方程求解下面设计一个求解线性方程组的函数文件line_solution.m。function[x,y]=line_solution(A,b)[m,n]=size(A);y=[];ifnorm(b)>0%非齐次方程组

ifrank(A)==rank([A,b])ifrank(A)==n%有唯一解

disp('原方程组有唯一解x')x=A\b;else%方程组有无穷多个解，基础解系

disp('原方程组有无穷个解，特解为x，其齐次方程组的基础解系为y')x=A\b;y=null(A,'r');

endelse

disp('方程组无解')%方程组无解

x=[];endelse%齐次方程组

disp('原方程组有零解x')x=zeros(n,1);%0解

ifrank(A)<ndisp('方程组有无穷个解，基础解系为y')%非0解

y=null(A,'r');endend第8章MATLAB方程求解例8-8求解方程组。程序如下：A=[1,-2,3,-1;3,-1,5,-3;2,1,2,-2];b=[1;2;3];[x,y]=line_solution(A,b)程序运行结果为：方程组无解x=[]y=[]说明该方程组无解。第8章MATLAB方程求解例8-9求方程组的通解。程序如下：formatrat%指定有理式格式输出A=[1,1,-3,-1;3,-1,-3,4;1,5,-9,-8];b=[1,4,0]';[x,y]=line_solution(A,b);x,yformatshort%恢复默认的短格式输出程序运行结果为：原方程组有无穷个解，特解为x，其齐次方程组的基础解系为y>Inline_solution(line11)Inaaa(line4)警告:秩不足，秩=2，tol=8.837264e-15。第8章MATLAB方程求解x=00-8/153/5y=3/2-3/43/27/41001所以原方程组的通解为：X=k1

+k2+，其中k1、k2为任意常数。第8章MATLAB方程求解8.2非线性方程数值求解非线性方程的求根方法很多，常用的有牛顿迭代法，但该方法需要求原方程的导数，而在实际运算中这一条件有时是不能满足的，所以又出现了弦截法、二分法等其他方法。在MATLAB中，非线性方程的求解和最优化问题往往需要调用最优化工具箱来解决。优化工具箱提供了一系列的优化算法函数，可用于解决工程中的最优化问题，包括非线性方程求解、极小值问题、最小二乘问题等。第8章MATLAB方程求解8.2.1单变量非线性方程求解在MATLAB中提供了一个fzero函数，可以用来求单变量非线性方程的根。该函数的调用格式为：z=fzero(filename,x0,tol,trace)其中，filename是待求根的函数文件名，x0为搜索的起点。一个函数可能有多个根，但fzero函数只给出离x0最近的那个根。tol控制结果的相对精度，默认时取tol=eps，trace指定迭代信息是否在运算中显示，为1时显示，为0时不显示，默认时取trace=0。第8章MATLAB方程求解例8-10求在x0=-5和x0=1作为迭代初值时的根。先建立函数文件fz.m。functionf=fz(x)f=x-1./x+5;然后调用fzero函数求根，命令如下：>>fzero(@fz,-5)%以-5作为迭代初值ans=-5.1926>>fzero(@fz,1)%以1作为迭代初值ans=0.1926可以绘制f(x)的图像（如下图所示），从中可以看出f(x)在x=-5和x=1附近的零点。第8章MATLAB方程求解第8章MATLAB方程求解8.2.2非线性方程组的求解在MATLAB的最优化工具箱中提供了非线性方程组的求解函数fsolve，其调用格式为：X=fsolve(filename,X0,option)其中，X为返回的解，filename是用于定义需求解的非线性方程组的函数文件名，X0是求解过程的初值，option用于设置优化工具箱的优化参数。优化工具箱提供了许多优化参数选项，用户在命令行窗口输入下列命令，可以将优化参数全部显示出来：>>optimset如果希望得到某个优化函数（例如fsolve函数）当前的默认参数值，则可在命令行窗口输入命令：>>optimsetfsolve如果想改变其中某个参数选项，则可以调用optimset函数来完成。例如，Display参数选项决定函数调用时中间结果的显示方式，其中'off'为不显示，'iter'表示每步都显示，'final'只显示最终结果。如果要将设定Display选项为off，可以使用命令：>>option=optimset('Display','off')。第8章MATLAB方程求解例8-11求下列方程组在（1，1，1）附近的解并对结果进行验证。首先建立函数文件myfun.m。functionF=myfun(X)x=X(1);y=X(2);z=X(3);F(1)=sin(x)+y+z^2*exp(x);F(2)=x+y+z;F(3)=x*y*z;在给定的初值x0=1，y0=1，z0=1下，调用fsolve函数求方程的根，命令如下：>>option=optimset('Display','off');>>X=fsolve(@myfun,[1,1,1],option)X=0.0224-0.0224-0.0000将求得的解代回原方程，可以检验结果是否正确，命令如下：>>q=myfun(X)q=1.0e-06*-0.5931-0.00000.0006第8章MATLAB方程求解例8-12求圆和直线的两个交点。圆：直线：该问题即为求解方程组：使用fsolve函数求解方程组时，必须先估计出方程组的根的大致范围。所给直线的方向数是(-12，24，-14)，故其与球心在坐标原点的球面的交点大致是(-1，1，-1)和(1，-1，1)。以这两点作为迭代初值。第8章MATLAB方程求解先建立方程组函数文件fxyz.m。functionF=fxyz(X)x=X(1);y=X(2);z=X(3);F(1)=x^2+y^2+z^2-9;F(2)=3*x+5*y+6*z;F(3)=x-3*y-6*z-1;再在MATLAB命令行窗口，输入如下命令：>>X1=fsolve('fxyz',[-1,1,-1],optimset('Display','off'))%求第一个交点X1=-0.95082.4016-1.5259>>X2=fsolve('fxyz',[1,-1,1],optimset('Display','off'))%求第二个交点X2=1.4180-2.33611.2377第8章MATLAB方程求解8.3最优化问题求解最优化方法包含有多个分支，如线性规划、整数规划、非线性规划、动态规划、多目标规划等。利用MATLAB的优化工具箱，可以求解线性规划、非线性规划和多目标规划问题。MATLAB还提供了求非线性函数最小值问题的求解方法，为优化方法在工程中的实际应用提供了更方便快捷的途径。第8章MATLAB方程求解8.3.1无约束最优化问题求解无约束最优化问题的一般描述为：其中，，该数学表示的含义亦即求取一组x，使得目标函数f(x)为最小，故这样的问题又称为最小化问题。MATLAB提供了3个求最小值的函数，它们的调用格式为：①[x,fval]=fminbnd(filename,x1,x2,option)：求一元函数在(xl，x2)区间中的极小值点x和最小值fval。②[x,fval]=fminsearch(filename,x0,option)：基于单纯形算法求多元函数的极小值点x和最小值fval。③[x,fval]=fminunc(filename,x0,option)：基于拟牛顿法求多元函数的极小值点x和最小值fval。MATLAB没有专门提供求函数最大值的函数，但只要注意到-f(x)在区间（a，b）上的最小值就是f(x)在（a，b）的最大值，所以fminbnd(-f,x1,x2)返回函数f(x)在区间（x1，x2）上的最大值。第8章MATLAB方程求解例8-13求函数在区间(-10，-1)和(1，10)上的最小值点。命令如下：>>f=@(x)x-1./x+5;>>[x,fmin]=fminbnd(f,-10,-1)%求函数在(-10，-1)内的最小值点和最小值x=-9.9999fmin=-4.8999>>fminbnd(f,1,10)%求函数在(1，10)内的最小值点ans=1.0001第8章MATLAB方程求解例8-14设求函数f在(0.5，0.5，0.5)附近的最小值。建立函数文件fxyz.m。functionf=fxyz0(u)x=u(1);y=u(2);z=u(3);f=x+y.^2./x/4+z.^2./y+2./z;在MALAB命令行窗口，输入如下命令：>>[U,fmin]=fminsearch(@fxyz0,[0.5,0.5,0.5])%求函数的最小值点和最小值U=0.50001.00001.0000fmin=4.0000第8章MATLAB方程求解8.3.2有约束最优化问题求解有约束最优化问题的一般描述为：其中，，该数学表示的含义亦即求取一组x，使得目标函数f(x)为最小，且满足约束条件G(x)≤0。记号s.t.是英文subjectto的缩写，表示x要满足后面的约束条件。约束条件可以进一步细化为:①线性不等式约束：Ax≤b。②线性等式约束：Aeqx=beq。③非线性不等式约束：C(x)≤0。④非线性等式约束：Ceq(x)=0。⑤x的下界和上界：Lbnd≤x≤Ubnd。第8章MATLAB方程求解MATLAB最优化工具箱提供了一个fmincon函数，专门用于求解各种约束下的最优化问题，其调用格式为：[x,fval]=fmincon(filename,x0,A,b,Aeq,beq,Lbnd,Ubnd,NonF,option)其中，x、fval、filename、x0和option的含义与求最小值函数相同。其余参数为约束条件，参数NonF为非线性约束函数的M文件名。如果某个约束不存在，则用空矩阵来表示。第8章MATLAB方程求解例8-15求解有约束最优化问题。首先编写目标函数M文件fop.m。functionf=fop(x)f=0.4*x(2)+x(1)^2+x(2)^2-x(1)*x(2)+1/30*x(1)^3;再设定约束条件，并调用fmincon函数求解此约束最优化问题，程序如下：x0=[0.5;0.5];A=[-1,-0.5;-0.5,-1];b=[-0.4;-0.5];lb=[0;0];option=optimset;option.LargeScale='off';option.Display='off';[x,f]=fmincon('fop',x0,A,b,[],[],lb,[],[],option)程序运行结果为：x=0.33960.3302f=0.2456第8章MATLAB方程求解8.3.3线性规划问题求解线性规划是研究线性约束条件下线性目标函数的极值问题的数学理论和方法。线性规划问题的标准形式为：在MATLAB中求解线性规划问题使用函数linprog，其调用格式为[x,fval]=linprog(f,A,b,Aeq,beq,Lbnd,Ubnd)其中，x是最优解，fval是目标函数的最优值。函数中的各项参数是线性规划问题标准形式中的对应项，x、b、beq、Lbnd、Ubnd是向量，A、Aeq为矩阵，f为目标函数系数向量。第8章MATLAB方程求解例8-16求解线性规划问题。

命令如下：f=[-5;-4;-6];A=[1,-1,1;3,2,4;3,2,0];b=[20;42;30];Aeq=[];Beq=[];LB=zeros(3,1);[x,favl]=linprog(f,A,b,Aeq,Beq,LB)程序运行结果为：Optimizationterminated.x=0.000015.00003.0000favl=-78.0000第8章MATLAB方程求解8.4常微分方程初值问题的数值求解考虑常微分方程的初值问题：y'=f(t,y)，t0≤t≤Ty(t0)=y0所谓其数值解法，就是求它的解y(t)在节点t0<t1<…<tm处的近似值y0，y1，…，ym的方法。所求得的y0，y1，…，ym称为常微分方程初值问题的数值解。一般采用等距节点tn=t0+nh，n=0，1，…，m，其中h为相邻两个节点间的距离，叫做步长。常微分方程初值问题的数值解法多种多样，比较常用的有欧拉（Euler）法、龙格—库塔（Runge-Kutta）法、线性多步法、预报校正法等。第8章MATLAB方程求解8.4.1龙格—库塔法简介对于一阶常微分方程的初值问题，在求解未知函数y时，y在t0点的值y(t0)=y0是已知的，并且根据高等数学中的中值定理，应有y(t0+h)=y1≈y0+hf(t0,y0)，h>0y(t0+2h)=y2≈y1+hf(t1,y1)一般地，在任意点ti=t0+ih，有：y(t0+ih)=yi≈yi-1+hf(ti-1,yi-1)，i=1，2，…，n当（t0，y0）确定后，根据上述递推式能计算出未知函数y在点ti=t0+ih，i=0，1，…，n的一列数值解：yi=y0，y1，y2，…，yn，i=0，1，…，n当然，递推过程中有一个误差累计的问题。在实际计算过程中，使用的递推公式一般进行过改造，著名的龙格—库塔公式是：y(t0+ih)=yi≈yi-1+(k1+2k2+2k3+k4)其中：k1=f(ti-1,yi-1)k2=f(ti-1+h/2,yi-1+h/2k1)k3=f(ti-1+h/2,yi-1+h/2k2)k4=f(ti-1+h,yi-1+hk3)第8章MATLAB