嗨~那么之前我们已经学习了（ahhh其实是我已经学习了）如何使用逐步逼近法来求解一个较麻烦但是满足利普希茨条件的初值问题的近似解。但是，现在情况是这样的，还记得之前我们的所求的函数是如何定义的吗？哈当然是定义在我们所给的矩形区域 
R：|x−x0|≤a,|y−y0|≤b
这个区域上面的了（如果你要研究的初值问题是没有给a和b的话，那么这个矩形区域R就需要靠你自己来定义了）。

那么既然我们的函数都定义在这个初值点的矩形区域内，那么自然，我们最后得到的结论，即该函数解的存在性和唯一性也一定是指的在这个矩形区域R内存在且唯一了。

那我们满足了吗？并没有，只在矩形区域R里面多憋屈呀，咱们开放一点，来一个在整个宇宙里都只有1个解的定理，多爽。嗯嗯，带着这样远大理想，我们来研究研究解的延拓。解的延拓是再说啥呢，就像我刚刚说的，你在这个小小的去建立存在唯一是不够好的，我要你在更大，更宽广的天空里存在且唯一。

为了更明确我们的目标，我们这里给出一个具体的问题供大家参考：


  
讨论方程dydx=y2−12的分别通过点(0,0),(ln2,−3)的解的存在区间。


这和以前的问题有什么不同，以前我们学会判断这个函数如果初值条件是这两个点的话，那么我们可以判断它在这两个点的附近的解是否存在，且唯一。但现在呢，好的，你用的你利普希茨条件判断出来了他们存在且唯一了，但是我现在要你说他们的解的存在区间是多少【注意，微分方程的解是一条积分曲线，所以这里其实问的是积分曲线的存在区间】

一脸懵逼有木有，而手握课本的我现在已经知道了该如何判断。但是，书上不给证明的阐述这条定理，哎，真不知道这样的读书有什么用，脑袋里装着一堆定理却缺少深邃的思维和自主学习的能力。完全是在凭借记忆力来获得好成绩。

我们这里当然会说明解的延拓的牛逼思想啦~

ok，那现在让我们开始，因为根据我之前所说，我们的起点也是在我们已经判断了某一只小函数它已经在初值位置的附近满足了存在唯一性的时候开始。

现在这里有一只小函数dydx=f(x,y)，它呢在一个小区间R里面满足了我们之前所说的利普希茨条件，那么如果这个小函数它还比较牛逼，不仅仅在R里面满足了利普希茨条件，而且它还在R的外面连续（就连续不断），那么它因为有辣么一点小牛逼，所以我们给它所具备的性能名称升个级，叫做局部的利普希茨条件【我怎么看这名字那么像降级了呢，它明明是变牛逼了的说。。。好吧让我们忽略这个问题】。


  
那么现在我们这个小函数就假设它已经满足了局部利普希茨条件了，记住它现在解的状态是在横坐标区间|x−x0|≤h的区间内存在且唯一，现在我们来考虑这个区间的边界时小函数f(x,y)的状态。这个小函数在边界上的这个点，也是满足利普希茨条件的（这是一个事实，因为他连续所以|f(x,y1)−f(x,y2)||y1−y2|这个值不可能达到无限大，所以就必然存在一个常数L使其满足利普希茨条件）【但是讲真，我对利普希茨条件中的y1和y2的理解还不够透彻呀~】。 

  但无论如何，一旦它边界上的点发生了满足利普希茨条件这种事的话，这意味着小函数的存在唯一性就可以扩展出去了呀！只要连续就能扩展，直至间断。


上述文字，说的大概就是延拓的思想吧。艾玛记定理有什么用，其延拓的思想才是有价值的东西呀，只关注定理定理无异于将思想取头取尾，而对其中的坎坷蜿蜒，对其中的峰回路转一概不看，就仿佛权利的游戏只看了开头和结局，就像人的一生只有出生和死亡一样。如果起点和终点是你在临走是所记得的东西的话，那该是多么的遗憾呀。

不扯了。。。我们还是来用定理总结一下刚刚我们所说，因为我特别讨厌定理这种表达方式，所以就不对它做强调处理了。

就和我们刚才所讨论的一样，如果有一个小函数它具有局部利普希茨条件这种较为牛逼的性质的话，那么介个小函数就可以通过在它连续范围G内任何与一个点的解开始延拓到无限接近于G的边界。

艾玛现在看来真是无语那些靠背定理，靠刷题来提高自己成绩，却对其中深邃的原理置之不理的人呀。与其像机器一样死记硬背，像工人一样通过重复提高熟练度，我真的觉得。哎，读什么大学，还不如去读技校呢···

继续的话，延拓的学习基本上已经说完了ahhh，在扯了一大堆淡之后，重要的只有那段原理吧ahhhh


  
补充说明： 

  emmm，说什么来着。。。就是只要连续就可以扩展直至间断点的出现是没错的，那间断点长什么样子呢？如果学过高等数学的同学应该对第一类间断点和第二类间断点还有印象吧。
  
  
简单来说第一类间断点就是某一个点的左右极限都存在，但是他们不相等，形象来说就是一条线被剪成两段， 躺在桌面上。 

  而第二类间断点就段得更彻底了，是这一点的极限不存在(±∞)
  
  
欸···为啥米要回顾第一类间断点呢，因为就算小函数它虽然本身看上去可能是在一个无界的范围内满足局部利普希茨条件，但是它的解y=φ(x)还是会受到第一间断点和第二间断点的约束，即一旦一条积分曲线在延拓的过程中遇到了第一或者第二类间断点，那么它的延拓生涯也就结束啦。就像我们接下来的实验一样。


解决实验

那我们现在就来搞定本次研究的实验吧~


  
讨论方程dydx=y2−12的分别通过点(0,0),(ln2,−3)的解的存在区间。
  
  
介个东西也没有给出a和b，但是我们不难发现这个f(x,y)是在某个区间上一定是连续且有界的。）所以它满足局部利普希茨条件，所以他满足在初值点附近的区间满足解的存在且唯一的性质，emmm，而且它满足局部利普希茨条件的区间是无界的，但是有我们刚才的补充说明，即使光看f的局部利普希茨条件是无界的，但是这个微分方程的解y=φ(x)也不一定能够在整个坐标系上都能够延拓，还是会受到第一类和第二类间断点的约束条件。所以接下来我们就按照书上的方法验证一下咯~
  
  
先求出这个方程的通解：
  
  
这个方程很显然是一个变量分离方程啦，也就是可以把x和y分离开，就像这样： 
1y2−1dy=12dx


  两边同时积分，左边就可以是： 
∫1y2−1dy=12∫1(y−1)−1(y+1)dy=12[ln(y−1)−ln(y+1)]+c1
右边就变成： 
∫12dx=x2+c2


  如果我们记c为任意常数。。。用c1,c2太痛苦了，辣么我们现在将左右两边相等，就有(轻微化简了一下)： 
lny−1y+1=x−c


  因为我们要求解的是y，所以我们继续化简 
y−1y+1=ex−c


  方程左边:y+1−2y+1=1−2y+1=右边=ex−c，化简得到： 
2y+1=1−ex−c


  所以有: 
y=21−ex−c−1=1+ex−c1−ex−c=1+c⋅ex1−c⋅ex


  这个就是它的通解了，然后如果要这个积分曲线经过(0,0)点的话，即c=−1，这时 
y=φ1(x)=1−ex1+ex
，显然肯定是没有第一类间断点的，而因为这个表达式的分母不可能等于0所以也不会有第二类间断点的。所以它的解的存在且唯一的区间就是(−∞,+∞)
那看完了它我们再来看看过(ln2,−3)的那个点的积分曲线的存在区间，将初值条件代入我们的通解y=1+c⋅ex1−c⋅ex，则可以解得c=1，所以也就是这条积分曲线是 
y=φ2(x)=1+ex1−ex


  这个样子的。还是不会有第一类间断点存在，当x从ln2向右延拓的时候，当然是一路无阻可以延拓到+∞，但是往左延拓的时候则会在x=0的时候出现第二类间断点，趋近于−∞，所以总结一下如果初值条件为(ln2,−3)的话，那么它的解的存在且唯一的区间就是(0,+∞)了。


emmm以上就是这个实验的完整版的解答了，如果有任何疑问或者改进的话，等有了再来继续补充吧~

算法原理和程序框图
龙格—库塔法是一种求其准确解y(x)在一系列点xi处y（xi）的近似值yi的方法，yi称为数值解。经典的四阶龙格库塔法方程如下：
y'=f(t,y),y(t0)=y0输出按如下求解yn+1=yn+h(k1+2k2+3k3+4k4)/6其中                       k1=f(tn,yn)
k2=f(tn+h/2,yn+hk1/2)
k3=f(tn+h/2,yn+hk2/2)
k4=f(tn+h,yn+hk3)
4.2程序使用说明
本程序单独编写了一阶初值问题单独求解和高阶常微分方程的求解问题。常微分方程的阶数，x的上下限和步长直接输入就好。常微分方程则在系统提示input后再输入，并且按照y(1,1)=y,y（2,1）=y’,y(3,1)=y’’的格式输入,然后输入初始条件f(x),f’(x),f’’(x)等计算结果。最后求得所要的函数值。
4.3计算实例展示
例题为课本304页计算实习9.2。计算y（1）的结果如下：
ans =
  -1.000000000000000
  -0.847436458333333
  -0.689482936121419
  -0.525724908067489
  -0.355719511320969
  -0.178993729355764
   0.004957535888930
   0.196673510288970
   0.396729719312157
   0.605740292781225
   0.824360410550626
   1.053288897488131
   1.293270976642268
   1.545101189994835
   1.809626496745704
   2.087749559656611
% ** 文件名:Runge_Kutta.m
% 
% ** 创建人:**
% 
% ** 日 期:2016.12.14
% 
% ** 描 述:四级四阶龙格－库塔法求解常微分方程的初值问题
% 
% ** 函数：第304 页计算实习 9.2   
% 
% ** 参考教材：《数值分析》李乃成，梅立泉，科学出版社
%
clear
clc
format long
% % % m=input('请输入常微分方程的阶数m=');
% % % a=input('请输入x下限a=');
% % % b=input('请输入x上限b=');
% % % h=input('请输入步长h=');
% % % ym=input('令y(1,1)=y,y（2,1）=y’,y(3,1)=y’’...请输入ym=','s');     %输入的时候必须按照这个形式输入y1=y(1,1);

m=3;
a=0;
b=1;
h=0.05;
%ym=y(3,1)+y(2,1)-y(1,1)+2*x-3;
ym=input('令y(1,1)=y,y（2,1）=y’,y(3,1)=y’’...请输入ym=','s');     %输入的时候必须按照这个形式输入y1=y(1,1);
% if m==1                                               %一阶初值问题单独求解
%     mm=(b-a)/h;
%     y(1,1)=input('请输入在初值点的函数值f(a)=');
%     x=a;
%     y11(1)=y(1,1);
%     for k1=2:(mm+1)
%         y1=y(1,1);
%         K(1,1)=h*(eval(ym));                         %计算K1
%         x=x+h/2;
%         y(1,1)=y1+K(1,1)/2;
%         y1=y(1,1);
%         K(1,2)=h*(eval(ym));                         %计算K2
%         x=x;
%         y(1,1)=y1+K(1,2)/2-K(1,1)/2;
%         y1=y(1,1);
%         K(1,3)=h*(eval(ym));                          %计算K3
%         x=x+h/2;
%         y(1,1)=y1+K(1,3)-K(1,2)/2;
%         y1=y(1,1);
%         K(1,4)=h*(eval(ym));                          %计算K4
%         y11(k1)=y11(k1-1)+(K(1,1)+2*K(1,2)+2*K(1,3)+K(1,4))/6;
%         y(1,1)=y11(k1);
%         x=a+(k1-1)*h;
%        
%     end
% y11
% else                                                 %高阶初值问题
    mm=(b-a)/h;                                     %一共要求解mm个数据点
    for k2=1:m                                      %读取初值条件
        fprintf('请输入%d阶导数的初值f(%d)(a)=\n',(k2-1),(k2-1));
        y(k2,1)=input('=');
    end
    for k2=1:m                                        
         y22(1,k2)=y(k2,1);                          %先把初值保存在矩阵y22(m,n)中,m表示第几个所求点，n表示第n阶初值
    end
%%%%初值

% y22(1,1)=-1;
% y22(1,2)=3;
% y22(1,3)=2;
    x=a;
    for k4=2:(mm+1)                                   %求解mm个数据点的循环
        for k=1:(m-1)                                 %计算K1，包括每一阶的K1
            K(k,1)=h*y(k+1,1);                        %y(k+1,1)中k+1表示第k+1阶，1表示第一个点；K(k,1)中k表示阶数，1表示K1
        end
        K(m,1)=h*(eval(ym));
        x=x+h/2;                                      %求解K1之前，先重新对x和y赋值
        for k3=1:m              
            y(k3,1)=y(k3,1)+K(k3,1)/2;
        end
        for k=1:(m-1)                                  %计算K2
            K(k,2)=h*y(k+1,1);
        end
        K(m,2)=h*(eval(ym));
        x=x;
        for k3=1:m
            y(k3,1)=y(k3,1)-K(k3,1)/2+K(k3,2)/2;
        end
        for k=1:(m-1)                                  %计算K3
            K(k,3)=h*y(k+1,1);
        end
        K(m,3)=h*(eval(ym));
        x=x+h/2;
        for k3=1:m
            y(k3,1)=y(k3,1)+K(k3,3)-K(k3,2)/2;  %这里容易出错
        end
        for k=1:(m-1)                                     %计算K4
            K(k,4)=h*y(k+1,1);
        end
        K(m,4)=h*(eval(ym));
        for k5=1:m
            y22(k4,k5)=y22(k4-1,k5)+(K(k5,1)+2*K(k5,2)+2*K(k5,3)+K(k5,4))/6;         %这里，除了要求出下一个点的数值，还要求出相应的导数值
        end
        for k6=1:m                                       %除了对y(1,1)重新赋值外，还要对y(2,1)等重新赋值
            y(k6,1)=y22(k4,k6); 
        end                                                      
        x=a+(k4-1)*h;
    end
        y22(:,1)
%end

工程上 ,不少实际问题的数学模型都涉及到非线性方程 f ( x )= 0的求解。由于工程问题对应的方程 f (x )= 0大多不存在求根公式 ,因此确定精确解十分困难。 故近似解的计算成为人们关心的主要问题。 目前 ,人们已提出了不少求解非线性方程 f (x )= 0近似解的方法 ,其中 ,牛顿迭代法是最基本的方法 之一。 由于牛顿迭代法在方程的单根附近具有平方收敛速度 ,而且还可以求解方程的重根和复根 ,故该法得到了广泛应用。但牛顿迭代法对初值的要求较苛刻 ,只有适当选取初值 ,才能保证其收敛性。

为保证牛顿迭代法局部收敛 ,必须对 f (x )和初值附加如下条件：

( 1) f (x )在区间 [a,b]上二阶可导; 

( 2) f (a) f (b) < 0; 

( 3) f′ (x )≠0; 

( 4) f″ (x )在区间 [a,b ]上不变号 ; 

( 5) 初值 x0∈ [a,b ]应使 f″ (x0 ) f ( x0 )> 0。

 其中 ,条件 ( 1)、 ( 2)保证了方程 f (x )= 0根的存在性;条件 ( 3)表示函数单调变化 ,则方程的根惟一;条件 ( 4)表示函数 f (x )的图形凸凹向不变;条件 ( 5)是对初值的要求

只要满足以上条件 ,就可保证牛顿迭代法的收敛性。但当表达式 f ( x )很复杂时 ,确定 f ( x )导数的计 算量很大 ,且求出 f″ (x )的表达式非常复杂 ,要满足条件 ( 3)、 ( 4)比较困难。另外 ,还有不少方程不完全符 合上述条件 ,无法按以上要求确定初值 ,迭代过程很难确保收敛性。

在文献一中提到了一种通过不断增加非线性方程对应曲线切线的斜率 ,达到迭代收敛的目的“改变斜率法”

 

在定义域内取任一初值 x0∈ [a,b ],按牛顿迭代法的公式进行计算 ,得



若 x1落在区间 [a,b]外 ,可将式 ( 1)中导数增加一倍 ,然后再重新进行迭代 ,得



若还落在区间 [a,b]外 ,可将导数再增加一倍 ,然后再重新计算 ,直到 x1 落在区间 [a,b ]内为止。 以下迭 代过程与此类似 ,其迭代格式为



参数 m的选取原则是使 xk+ 1∈ (a,b)。以上思路的几何意义如图 1所示。当曲线 y= f (x )在点 xk 处 的切线 (图 1中的实线 )与坐标轴的交点 (该点的横坐标即为 xk+ 1 )落在区间 [a,b]外时 ,将切线斜率增加 一倍 ,即使其变陡 ,成为割线 ,如图 1中的虚线所示。 这样 ,其与坐标轴的交点将落在区间 [a,b ]内 ,从而 保证迭代收敛。若该割线与坐标轴的交点仍然落在区间 [a,b]外 ,则使其斜率再增加一倍 ,直到满足要求 为止。 对于图 2所示的情形 ,显然不满足上述牛顿迭代法的局部收敛条件 ( 4) (因为 f″ ( x )正负交替变号 ) , 很难选取能保证收敛的合适初值。但采用本文提出的改变斜率法 ,取区间 [a,b ]内的任意初值 ,均可保证 迭代收敛。



图 1　改变斜率法的几何意义 图 2　异形曲线的迭代 　　需要指出的是 ,在迭代过程中 ,每一步迭代都要检查迭代的根是否落在区间 [a,b]外 ,若是 ,则将斜 率增加一倍重新计算即可。 显然 ,m= 0即为常规牛顿迭代法 ,所以常规牛顿迭代法可视为本文方法的特例。

 

 
试着用牛顿迭代法求 的根




clear,clc
x=0.1;%zm初始值
xs=[x];
n=30;%迭代次数
for i=1:n
    dzm=3*x;
    x=x-dzm;
    xs=[xs,x];
end
plot(xs,'*')

当m=0时，迭代所得的方程的根如下图



可以看到迭代结果是发散的

但是当m=1时，迭代结果如下图，为收敛的结果



这时的结果则是我们所需要的。

试求方程 f (x )= 5x^3- x^2- 1= 0在区间 [0. 2, 1. 5]内的根

clear,clc
x=0.1;%zm初始值
xs=[x];
result=0;
n=30;%迭代次数
for i=1:n
    dzm=(5*x^3-x^2-1)/(15*x^2-2*x);
    x=x-dzm;
    xs=[xs,x];
    a=xs(length(xs)-1);
    b=xs(length(xs));
    rab=abs((b-a)/b);
    if rab<1e-4
        result=b;
        break
    end
end
disp(result)

用牛顿迭代法迭代出的结果如下图，所得结果收敛





[1]陈玉骥.牛顿迭代法的一种改进方法[J].佛山科学技术学院学报(自然科学版),2012,30(05):1-3.

[2] 牛顿迭代的详细原理，以及初值对迭代结果的影响可以参照https://blog.csdn.net/ccnt_2012/article/details/81837154

一、基本概念

1. 数值解的概念

数值解(numerical solution)是指在特定条件下通过近似计算得出来的一个数值，是采用某种计算方法，如有限元的方法、 数值逼近、插值的方法 得到的解。

例如，求微分方程的初值问题

                                                                      

的解，可以从初值条件出发，按照一定的步长h，依据某种方法逐步计算微分方程解y(x)的近似解，这里 。这样求出的解称为数值解。

2. 解析解的概念

解析解(analytical solution)，就是给出解的具体函数形式，从解的表达式中就可以算出任何对应值。即包含分式、三角函数、指数、对数甚至无限级数等基本函数的解的形式。给出解的具体函数形式，从解的表达式中就可以算出任何对应值。用来求得解析解的方法称为解析法，解析法是常见的微积分技巧，如分离变量法等。

 

二、数值解法优势及常用方法

1. 欧拉法

 

2. 龙格-库塔法

龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。

该算法是构建在泰勒级数的基础之上的。对于一阶精度的拉格朗日中值定理有：

对于微分方程：y'=f(x,y)

y(i+1)=y(i)+h*K1

           K1=f(xi,yi)

当用点xi处的斜率近似值K1与右端点xi+1处的斜率K2的算术平均值作为平均斜率K*的近似值，那么就会得到二阶精度的改进拉格朗日中值定理：

y(i+1)=y(i)+[h*( K1+ K2)/2]

           K1=f(xi,yi)

           K2=f(x(i)+h,y(i)+h*K1)

依次类推，如果在区间[xi,xi+1]内多预估几个点上的斜率值K1、K2、……Km，并用他们的加权平均数作为平均斜率K*的近似值，显然能构造出具有很高精度的高阶计算公式。经数学推导、求解，可以得出四阶龙格－库塔公式，也就是在工程中应用广泛的经典龙格－库塔算法：

y(i+1)=y(i)+h*( K1+ 2*K2 +2*K3+ K4)/6

           K1=f(x(i),y(i))

           K2=f(x(i)+h/2,y(i)+h*K1/2)

           K3=f(x(i)+h/2,y(i)+h*K2/2)

           K4=f(x(i)+h,y(i)+h*K3)

通常所说的龙格-库塔法是指四阶而言的，我们可以仿二阶、三阶的情形推导出常用的标准四阶龙格-库塔法公式。

三、解析解法优势及常用方法