牛顿下山法求非线性方程组解
牛顿法、简化牛顿法与牛顿下山法、弦截法、解非线性方程组的牛顿法
河北联合大学第7章 非线性方程组的数值解法§7.3 牛顿法 §7.4 简化牛顿法与牛顿下山法§7.5 弦截法 §7.6 解非线性方程组的牛顿法1. 什么是求解f x =0的牛顿法？它是否总是收敛的？若f *x =0，x *是单根，f 光滑，证明牛顿法是局部二阶收敛的。
答：按式x 1 n =x n —n n x f x f '(n=0,1,2……n )求方程f x =0的近似解的方法称为牛顿法；牛顿法不总是收敛的，它是局部收敛的；设函数()f x 在其零点*x 邻近二阶连续可微，且*()0f x ﾢﾢ，则存在0d >，使得对任意**0[,]x x x d d - �，Newton 法所产生的序列{}n x 至少二阶收敛于*x 。
证明 由1() (0,1,2,)()n n n n f x x x n f x =-=ﾢL 知迭代函数为()()()f x x x f x j =-ﾢ，且有2()()()[()]f x f x x f x j ﾢﾢﾢ=ﾢ，若()f x ﾢﾢ在*x 邻近连续，则()x j ﾢ在*x 邻近连续，且****2()()()01[()]f x f x x f x j ﾢﾢﾢ==<ﾢ当迭代函数()x j在*x 邻近有r 阶连续导数，且**=()x x j ,()*()0k x j =(1,,1)k r =-L ，0)(*)( x r j 则迭代序列{}n x 在点*x 邻近是r 阶收敛的。
可知Newton 法产生的迭代序列{}n x 至少二阶收敛于*x 。
2. 什么是弦截法？试从收敛阶及每步迭代计算量与牛顿法比较其差别。
答：弦截法是函数逼近法的一种，基本思想是用区间 x x kk ,1-上的割线近似代替目标函数的导函数的曲线。
并用割线与横轴交点的横坐标作为方程根的近似。
在Newton 迭代公式中，每次计算导数运算量很大，为了避免计算导数值,用差商代替导数)(x k f,得到迭代公式 按如下迭代公式计算方程的近似解称为弦截法。
Newton迭代法求解非线性方程
Newton迭代法求解非线性方程一、 Newton 迭代法概述构造迭代函数的一条重要途径是用近似方程来代替原方程去求根。
因此，如果能将非线性方程f (x )=0用线性方程去代替，那么，求近似根问题就很容易解决，而且十分方便。
牛顿(Newton)法就是一种将非线性方程线化的一种方法。
设k x 是方程f (x )=0的一个近似根，把如果)(x f 在k x 处作一阶Taylor 展开，即：)x x )(x ('f )x (f )x (f k k k -+≈ (1-1)于是我们得到如下近似方程：0)x x )(x ('f )x (f k k k =-+ (1-2)设0)('≠k x f ，则方程的解为：x ̅=x k +f (x k )f (x k )́(1-3)取x ~作为原方程的新近似根1+k x ，即令： )x ('f )x (f x x k k k 1k -=+, k=0,1,2,…(1-4)上式称为牛顿迭代格式。
用牛顿迭代格式求方程的根的方法就称为牛顿迭代法，简称牛顿法。
牛顿法具有明显的几何意义。
方程：)x x )(x ('f )x (f y k k k -+= (1-5)是曲线)x (f y =上点))x (f ,x (k k 处的切线方程。
迭代格式(1-4)就是用切线式(1-5)的零点来代替曲线的零点。
正因为如此，牛顿法也称为切线法。
牛顿迭代法对单根至少是二阶局部收敛的，而对于重根是一阶局部收敛的。
一般来说，牛顿法对初值0x 的要求较高，初值足够靠近*x 时才能保证收敛。
若要保证初值在较大范围内收敛，则需对)x (f 加一些条件。
如果所加的条件不满足，而导致牛顿法不收敛时，则需对牛顿法作一些改时，即可以采用下面的迭代格式：)x ('f )x (f x x k k k 1k λ-=+,⋯=,2,1,0k (1-6)上式中，10<λ<，称为下山因子。
牛顿下山法求非线性方程组解
《MATLAB 程序设计实践》课程考核1、编程实现以下科学计算算法，并举一例应用之。
（参考书籍《精通ＭＡＬＡＢ科学计算》，王正林等著，电子工业出版社，２００９年）“牛顿下山法求非线性方程组解”解：算法说明：牛顿下山法的迭代公式：())()(11n n n n x F x F w x x -+-=w 的取值范围为10≤w ，为了保证收敛，还要求w 的取值使得：())(1n n x F x F +可以用作次减半法来确定w 。
为了减少计算量，还可以用差商来代替偏导数。
在MA TLAB 中编程实现的非线形方程组的牛顿下山法的函数：mulDNewton 。
功能：用牛顿下山法求非线形方程组的一个解。
调用格式：[]),0(,eps x mulDNewton n r = 其中， x0为初始迭代向量Eps 为迭代精度；r 为求出的解向量n 为迭代步数。
牛顿下山法的MATLAB 代码如下：function [r,n]=mulDNewton(x0,eps) %牛顿下山法求非线形方程组的一个解 %初始迭代向量：x0 %迭代精度：eps %解向量： r %迭代步数：n if nargin==1 eps=1.0e-4;%输入的自变量的数目为1个时，精度定为eps=1.0e-4 endr=x0-myf(x0)/dmyf(x0); %当n=1时，取w=1 n=1; tol=1;%初始n 和tol 的值 while tol>eps x0=r; ttol=1; %初始ttol 的值 w=1;%初始w 的值，w 就是下山因子alpha F1=norm(myf(x0)); while ttol>=0r=x0-w*myf(x0)/dmyf(x0); ttol=norm(myf(r))-F1; w=w/2; endtol=norm(r-x0); n=n+1; if (n>100000)disp('迭代步数太多，可能不收敛！'); return ; end end举例说明：牛顿下山法求下面方程组⎩⎨⎧=--=-+0sin 1.0cos 5.00)cos(1.0sin 5.02211211x x x x x x x 的根，其初始迭代值取（0,0）。
解非线性方程组的牛顿迭代法
为克服这两个缺点，通常可用下述方法.
（1） 简化牛顿法，也称平行弦法.
xk 1 xk Cf ( xk )
其迭代公式为 （4.7）
C 0,1 ,.
迭代函数 ( x) x Cf ( x).
若在根 x * 附近成立 ( x) 1 Cf ( x) 1 ，即取 0 Cf ( x) 2，则迭代法（4.7）局部收敛.
8
xk
C 2 C
q2
k
1 q
2k
.
对任意 x0 0，总有 q 1，故由上式推知，当 k 时 xk C ，即迭代过程恒收敛. 例8 解 求 115 .
表7 6 计算结果 k 0 1 2 3 4 xk 10 10.750000 10.723837 10.723805 10.723805
f ( x) , f ( x)
由于
( x)
f ( x) f ( x) . 2 [ f ( x)]
假定 x *是 f ( x) 的一个单根，即 f ( x*) 0, f ( x*) 0 ， 则由上式知 ( x*) 0 ，于是依据定理4可以断定，牛顿法 在根 x *的邻近是平方收敛的.
准备 迭代
x0 ，计算 f 0 f ( x0 ), 选定初始近似值
步骤2
按公式
x1 x0 f 0 / f 0
迭代一次，得新的近似值 x1，计算 f1 f ( x1 ), f1 f ( x1 ). 步骤3 控制
x1 满足 1 如果
f1 2 ，则终 或
5
止迭代，以 x1作为所求的根；否则转步骤4. 允许误差，而
3
又因
( x*)
f ( x*) , f ( x*)
非线性方程的数值解法牛顿下山法matlab
非线性方程的数值解法——计算物理实验作业九陈万 物理学2013级 130********● 题目：用下列方法求0133=--=x x f(x)在20=x 附近的根。
根的准确值 87938524.1*=x ,要求计算结果精确到四位有效数字。
（1）用牛顿法；（2）用弦截法，取；9.1,210==x x● 主程序：clearclc;%----------------初值设定-------------------x0 = 2;x1 = 1.9;eps = 0.00001;N = 50;%----------------迭代求解-------------------Newton(x0,eps,N);Newton_downhill(x0,eps,N);Secant_Method(x0,x1,eps,N);● 子程序：f(x)function [y]=f(x)y = x^3-3*x-1; %函数f(x)End● 程序一：牛顿法function Newton(x0,eps,N)% 牛顿法% x0是迭代初值，eps是精度，N是迭代上限format long;k = 1;while(1)ff = (f(x0+0.1*eps)-f(x0))/(0.1*eps);if ff == 0disp('分母为零，请重新选择初始迭代值')break;elsex1=x0-f(x0)/ff ;if abs(x1-x0)<epsdisp('满足精度要求的根是：')disp(x1)break;elseif k>=Ndisp('迭代失败，请检查程序是否有误')breakelsek = k+1;x0 = x1;endendend程序二：弦截法function Secant_Method(x0,x1,eps,N)% 弦截法% x0,x1是迭代初值，eps是精度，N是迭代上限format long;k = 1;while(1)if f(x0)==0disp('满足精度要求的解是：')disp(x0)elseif f(x1)==0disp('满足精度要求的解是：')disp(x1)break;elseif abs(f(x1)-f(x0))==0disp('分母为零，请重新选择初始迭代值')break;elsex2 = x1-f(x1)*(x1-x0)/(f(x1)-f(x0));if abs(x2-x1)<epsdisp('满足精度要求的解是：')disp(x2)break;elseif k>=Ndisp('迭代失败，请检查程序是否有误')break;elsek = k+1;x0 = x1;x1 = x2;endendend程序三：牛顿下山法function Newton_downhill(x0,eps,N)% 牛顿下山法% x0是迭代初值，eps是精度，N是迭代上限format long;k = 1;while(1)lamda = 1;ff = (f(x0+0.1*eps)-f(x0))/(0.1*eps);if ff == 0disp('分母为零，请重新选择初始迭代值')elsewhile(1)x1 = x0-lamda*f(x0)/ff ;if f(x1)>=f(x0)lamda = 0.5*lamda;elsebreak;endendif abs(x1-x0)<epsdisp('牛顿下山法满足精度要求的根是：')disp(x1)break;elseif k>=Ndisp('迭代失败，请检查程序是否有误')breakelsek = k+1;x0 = x1;endendendend●程序运行结果：牛顿法：满足精度要求的根是：1.879385241571819弦截法：满足精度要求的解是：1.879385241572444●分析讨论：从运行结果来看，牛顿法与弦截法的结果与给定准确值完全相等；从运行时间上看速度都相当快。
牛顿拉夫逊法
并将雅可比矩阵分块，把每个2 2 阶子
阵
如MH
ij ij
Nij H ij
Lij
Rij
N ij S ij
等
作为分块矩阵的
元素，则按节点号顺序而构成的分块雅
可比矩阵将和节点导纳矩阵具有同样的
稀疏结构，是一个高度稀疏的矩阵。
(4) 和节点导纳矩阵具有相同稀疏结 构的分块雅可比矩阵在位置上对称，但 由于H ij H ji , Nij N ji , M ij M ji , Lij L ji ， 所以雅可比矩阵不是对称阵。
由上式根据初值 x(0) 可求得第一次迭 代的修正量
x(0) [ f (x(0) )]1 f (x(0) ) (1-25)
将 x(0) 和 x(0) 相加，得到变量的第一次 改进值 x(1) 。
因此，应用牛顿法求解的迭代格式为
f (x(k ) )x(k ) f (x(k) )
(1-26)
x(k1) x(k ) x(k )
（2）修正xi(1)＝ xi(0)＋ △xi(0) ，算出△f，J中各元素，代入 上式方程组，解出△ xi(1) ；
直至 f (x(k) ) 1 或 x(k) 2
注意：xi的初值要选得接近其精确值，否则将不迭代。
及 2n 1个，在 PV 节点比例不大时，两 者的方程式数目基本接近 2n 1个。
(2) 雅可比矩阵的元素都是节点电压 的函数，每次迭代，雅可比矩阵都需要 重新形成。
(3) 从雅可比阵非对角元素的表示式
可见，某个非对角元素是否为零决定于
相应的节点导纳矩阵元素 Yij 是否为零。
如将修正方程式按节点号的次序排列，
线性 迭代公式：
xk
1
xk
非线性方程(组)的数值解法——牛顿法、弦切法
(3) 用 Newton 法解 (x) = 0
x ( x 2 2) 3 ( x) x x2 2
ex76.m
14
弦截法与抛物线法
弦截法与抛物线法
目的：避免计算 Newton 法中的导数，且具有较 高的收敛性（超线性收敛） 弦截法（割线法）：用差商代替微商 抛物线法：用二次多项式近似 f(x)
2
x
k
C
2
2
xk 1 C xk C xk 1 C xk C 2k xk C x0 C xk C x0 C k q2 xk C 2 C 2k 1 q
q
2k
对任意 x0>0， 总有 |q|<1， 即牛顿法收敛
8
牛顿法
牛顿的优点
至少二阶局部收敛，收敛速度较快，特别是当迭代点 充分靠近精确解时。
牛顿法是目前求解非线性方程 (组) 的主要方法 牛顿的缺点
对重根收敛Βιβλιοθήκη 度较慢（线性收敛） 对初值的选取很敏感，要求初值相当接近真解 先用其它算法获取一个近似解，然后使用牛顿法
需要求导数！
9
简化的Newton法
f ( xk ) f '( xk ) 迭代格式： xk 1 xk [ f '( xk )]2 f ( xk ) f ''( xk )
13
举例
例：求 x4 - 4x2 ＋ 4=0 的二重根 x* 2 (1) 普通 Newton 法
x2 2 1 ( x ) x 4x
(2) 改进的 Newton 法 x2 2 2 ( x) x
简化的 Newton 法
135-7-4牛顿下山法、弦截法、解非线性方程组的牛顿法
4. -1.
计
Clear[x]
5. 0.
FindRoot[f[x]==0,{x,0.1I}]
-0.170616
{x -> -0.5 + 0.866025 I}
f ( x) x2 x 1 的图形
程 序 设 计
典型例题
典型例题
例1
用弦截法解方程
f ( x) x 3 7.7x 2 19.2x 15.3的根，取x0 1.5, x1 4
为了叙述方便，以解二阶非线性方程组为例演示解题方法和步骤
非 线
设二
阶
方
程
组：gf((xx,,
y) y)
0 0
性 方 程
写
成向量
形
式 ：F
(w)
f (x, g( x,
yy)),
其
中w
x y
组 的 牛
将f ( x, y)，g( x, y)在( x0 , y0 )附近作二元泰勒展开， 并取其线性部分，得到方程组：
x
f
( x1 , x
y1 )
y
f
(x1 , y
y1 )
f
( x1 ,
y1 )
x
g( x1 , x
y1
)
y
g( x1 , y
y1
)
g(
x1 ,
y1
)
方
f f
程 组 的 牛
记
J
x g
x
y g
y ( x1 , y1 )
当 J 0时，解出x, y
顿 解 法
迭 代：
w2
w1
x y
记作
x2 y2
依 此 类 推 可 得 ： w k 1
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
《MATLAB 程序设计实践》课程考核1、编程实现以下科学计算算法，并举一例应用之。
（参考书籍《精通ＭＡＬＡＢ科学计算》，王正林等著，电子工业出版社，２００９年）“牛顿下山法求非线性方程组解”解：算法说明：牛顿下山法的迭代公式：())()(11n n n n x F x F w x x -+-=w 的取值范围为10≤w ，为了保证收敛，还要求w 的取值使得：())(1n n x F x F +可以用作次减半法来确定w 。
为了减少计算量，还可以用差商来代替偏导数。
在MA TLAB 中编程实现的非线形方程组的牛顿下山法的函数：mulDNewton 。
功能：用牛顿下山法求非线形方程组的一个解。
调用格式：[]),0(,eps x mulDNewton n r = 其中， x0为初始迭代向量Eps 为迭代精度；r 为求出的解向量n 为迭代步数。
牛顿下山法的MATLAB 代码如下：function [r,n]=mulDNewton(x0,eps) %牛顿下山法求非线形方程组的一个解 %初始迭代向量：x0 %迭代精度：eps %解向量： r %迭代步数：n if nargin==1 eps=1.0e-4;%输入的自变量的数目为1个时，精度定为eps=1.0e-4 endr=x0-myf(x0)/dmyf(x0); %当n=1时，取w=1 n=1; tol=1;%初始n 和tol 的值 while tol>eps x0=r; ttol=1; %初始ttol 的值 w=1;%初始w 的值，w 就是下山因子alpha F1=norm(myf(x0)); while ttol>=0r=x0-w*myf(x0)/dmyf(x0); ttol=norm(myf(r))-F1; w=w/2; endtol=norm(r-x0); n=n+1; if (n>100000)disp('迭代步数太多，可能不收敛！'); return ; end end举例说明：牛顿下山法求下面方程组⎩⎨⎧=--=-+0sin 1.0cos 5.00)cos(1.0sin 5.02211211x x x x x x x 的根，其初始迭代值取（0,0）。
解：首先建立myf.m 函数文件，输入以下内容：function f=myf(x)f(1)=0.5*sin(x(1))+0.1*cos(x(2)*x(1))-x(1); f(2)=0.5*cos(x(1))-0.1*sin((2))-x(2); f=[f(1) f(2)];再建立dmyf.m 导数的雅克比矩阵，输入以下内容：function df=dmyf(x)df=[0.5*cos(x(1))-0.1*x(2)*sin(x(2)*x(1))-1,-0.1*x(1)*sin(x(2)*x(1));-0.5*sin(x(1)),-0.1*cos(x(2))-1];然后，再MATLAB 命令窗口中输入：[r,n]=mulDNewton([0,0]) 输出计算结果为： r =0.1981 0.3993n =5由计算结果，初始迭代值取（0，0）时，用5步迭代得到了方程组的一组解（0.1981，0.3993）。
运行结果如下图所示：流程图2．编程解决以下科学计算和工程实际问题。
1)在平炉炼钢过程中，由于矿石及炉气的氧化作用，铁水中的总含碳量在不断降低，一炉钢的冶炼初期（熔化期）中总的去碳量y于所加的两种矿石（天然矿石和烧结矿石）的量x1，x2及熔化时间x3有关（熔化时间愈长则去碳量愈多）。
经实测，某号平炉的49炉数据如下（在此省略），试求出y与x1，x2，x3之间的线性回归方程。
算法说明：在MA TLAB统计工具箱中使用函数regress实现多元线性回归。
功能：1确定几个特定的变量之间是否存在相关关系，如果存在的话，找出它们之间合适的数学表达式；2根据一个或几个变量的值，预测或控制另一个变量的取值，并且可以知道这种预测或控制能达到什么样的精确度；3进行因素分析。
例如在对于共同影响一个变量的许多变量（因素）之间，找出哪些是重要因素，哪些是次要因素，这些因素之间又有什么关系等等。
调用格式：b=regress(y,x)[b,bint,r,rint,stats]=regress(y,x,alpha)解：在MATLAB中编写程序dyxxhg.m来实现，具体代码如下：%dyxxhg.m 多元线形回归%%%%%输入数据x1=[2,7,5,12,1,3,2,6,7,0,3,0,8,6,0,3,7,16,6,0,9,4,0,9,2,9,12,6,12,0,5,4,0,6,4,10,4, 5,9,6,5,5,8,2,7,4,10,3,4];x2=[18,9,14,3,20,12,17,5,8,23,16,18,4,14,21,14,12,0,16,15,0,6,17,0,16,6,5,13,7,24,1 2,15,20,16,17,4,14,13,8,13,8,11,6,13,8,10,5,17,15];x3=[50,40,46,43,64,40,64,39,37,55,60,49,50,51,51,51,56,48,46,52,40,32,47,44,39,39,6 1,41,47,61,37,48,45,42,48,48,36,36,51,54,100,44,63,55,50,45,40,64,72];y=[4.3302,3.6485,4.4830,5.5468,5.4970,3.1125,5.1182,3.8759,4.6700,4.9536,5.0060,5.2 701,5.8772,5.4849,4.5960,5.6645,6.0795,3.2194,5.8076,4.7306,4.6805,3.1272,2.6104,3. 7174,3.8946,2.7066,5.6314,5.8152,5.1302,5.3910,4.4533,4.6569,4.5212,4.8650,5.3566,4 .6096,2.3815,3.8746,4.5919,5.1588,5.4373,3.9960,4.3970,4.0622,2.2905,4.7115,4.5310, 5.3637,6.0771];%%%%%%保存数据save dyxxhg_data x1x2x3y%%%%取出数据load dyxxhg_data%%%%执行回归命令x=[ones(49,1),x1',x2',x3'];[b,bint,r,rint,stats]=regress(y',x);b %输出参数bint %各参数估计的置信区间 stats %几个特殊统计量执行后，输出如下结果： >> dyxxhg b =0.7753 0.1504 0.1014 0.0370bint =-0.9819 2.5325 0.0270 0.2739 0.0245 0.1783 0.0157 0.0583stats =0.3331 7.4925 0.0004 0.6809所以，回归模行方程是：x x x y3210370.01014.01504.07753.0ˆ+++=此外，由stats 的值可知：R 2＝0.3331,F =7.4925,P =0.6809 运行结果如下图所示 ：流程图开始输入各变量的数值保存，加载数据执行回归命令，使用函数regress输出需求的相关参数 b ，bint，stats结束2) 用四阶Runge-Kutta 方法求解初值问题：101/10y(0)22≤≤⎪⎩⎪⎨⎧=-='-t yy e t取h ＝0.2，h ＝0.1，并与精确解y （0.4），y(1)比较。
其中精确解的解析表达式为：e et t t y 22101--+=。
算法说明：在龙格-库塔法中，四阶龙格-库塔法的局部截断误差约为o(h5)，被广泛应用于解微分方程的初值问题。
其算法公式为：)22(63211k k k hy y n n +++=+ 其中：⎪⎪⎪⎩⎪⎪⎪⎨⎧++=++=++==) ,()21,21()21 ,21(),(3423121hk y h x f k hk y h x f k hk y h x f k y x f k n n n n n n n n 功能：用四阶龙格-库塔算法求解常微分方程。
调用格式：[x,y]=MyRunge_Kutta(fun,x0,xt,y0,PointNum,varargin)。
其中，fun 为函数名； x0为自变量区间初值； xt 为自变量区间终值； y0为函数在x0处的值；PointNum 为自变量在[x0,xt]上所取的点数； Varargin 为可选项参数；x 为输出自变量区间上所取点的x 值； y 为对应点上的函数值。
解：原代码function [x,y] = MyRunge_Kutta(fun,x0,xt,y0,PointNum,varargin) %Runge-Kutta 方法解微分方程形为 y'(t)=f(x,y(x))%此程序可解高阶的微分方程。
只要将其形式写为上述微分方程的向量形式 %函数 f(x,y): fun%自变量的初值和终值：x0, xt%y0表示函数在x0处的值，输入初值为列向量形式 %自变量在[x0,xt]上取的点数：PointNum %varargin 为可输入项，可传适当参数给函数f(x,y) %x ：所取的点的x 值%y：对应点上的函数值if nargin<4 | PointNum<=0PointNum=100;endif nargin<3y0=0;endy(1,:)=y0(:)'; %初值存为行向量形式h=(xt-x0)/(PointNum-1); %计算步长x=x0+[0:(PointNum-1)]'*h; %得x向量值for k=1:(PointNum) %迭代计算f1=h*feval(fun,x(k),y(k,:),varargin{:});f1=f1(:)'; %得公式k1f2=h*feval(fun,x(k)+h/2,y(k,:)+f1/2,varargin{:});f2=f2(:)'; %得公式k2f3=h*feval(fun,x(k)+h/2,y(k,:)+f2/2,varargin{:});f3=f3(:)'; %得公式k3f4=h*feval(fun,x(k)+h,y(k,:)+f3,varargin{:});f4=f4(:)'; %得公式k4y(k+1,:)=y(k,:)+(f1+2*(f2+f3)+f4)/6; %得y(n+1)end在MA TLAB中编写程序rk.m来实现，具体代码如下：%运行四阶R-K法clear, clc %清除内存中的变量x0=0;xt=1;Num=100;h=(xt-x0)/(Num-1);x=x0+[0:Num]*h;fun =inline('exp(-2*x)-2*y','x','y'); %用inline构造函数f(x,y)y0=1/10; %设定函数初值PointNum=5; %设定取点数[x1,y1]=ode23(fun,[0,1],1/10);[xr,yr]=MyRunge_Kutta(fun,x0,xt,y0,PointNum);MyRunge_Kutta_x=xr'MyRunge_Kutta_y=yr'执行后，输出如下结果：MyRunge_Kutta_x =0 0.2500 0.5000 0.7500 1.0000MyRunge_Kutta_y =0.1000 0.2121 0.2205 0.1895 0.1488 0.1108>>y(0.4)=0.224664482,y(1)=0.148868811，相比较2者有所差距。