求几个函数，看课件上，老师推导下面这个函数，注意f(tn,phin)=dphi(tn)/dt(n)，乘一个M之后，Mdphi(tn)/dt(n)=B(tn)phi(tn)+C(tn), B和C的值可以参考老师课件上给的案例，还有要注意给的B的值，是[B(tn),B(tn+1)],还是[B(tn+1),B(tn)]
 
function [LHS,RHS] = OST(theta,timestep,M,B,C,sol)
LHS=M-theta*timestep*B(1);%[B(tn+1),B(tn)]
RHS=(M+(1-theta)*timestep*B(2))*sol+timestep*(theta*C(1)+(1-theta)*C(2));
end
 
根据老师推导的函数方法，同理可以推出其他的函数表达式
 
function [LHS,RHS] = AB2(timestep,M,B,C,sol)
LHS=M;
RHS=(M+3/2*timestep*B(1))*sol(1)+timestep/2*(3*C(1)-B(2)*sol(2)-C(2));%[B(tn),B(tn-1)],
end
 
function [LHS,RHS] = AM3(timestep,M,B,C,sol)
LHS=M-5/12*timestep*B(1);
RHS=(M+2/3*timestep*B(2))*sol(1)+timestep/12*(5*C(1)+8*C(2)-B(3)*sol(2)-C(3));
end
 
function [LHS,RHS] = BDF2(timestep,M,B,C,sol)
LHS=3/2*M-timestep*B;
RHS=2*M*sol(1)-1/2*M*sol(2)+timestep*C;
end
 
clear;
te=0:0.01:2;t=0:0.1:2;p1(1)=0;p2(1)=0;p3(1)=0;dt=0.1;p=exp(-5*te).*(te.*te-2*te+2)-2*exp(-6*te);M=1;C=t.*t.*exp(-5*t);%te和p是原函数的值

for i=1:length(t)-1%要求出所有的值，到i+1就是所有的值
    p1(i+1)=p1(i)+dt*(t(i)*t(i)*exp(-5*t(i))-6*p1(i));
    p2(i+1)=(p2(i)+dt*t(i+1)*t(i+1)*exp(-5*t(i+1)))/(1+6*dt);%注意这边要得到化简，i+1放一边
    p3(i+1)=(p3(i)+1/2*dt*(t(i+1)*t(i+1)*exp(-5*t(i+1))+t(i)*t(i)*exp(-5*t(i))-6*p3(i)))/(1+1/2*dt*6);%这是三种直接求解初值问题的方法
end
sol1=[p3(2) p3(1)];sol2=[p3(2) p3(1)];sol3=[p3(2) p3(1)];p4=[p3(1) p3(2)];p5=[p3(1) p3(2)];p6=[p3(1) p3(2)];B1=[-6 -6];B2=[-6 -6 -6];B3=[-6];
for j=2:length(t)-1
    [L1,R1]=AB2(dt,1,B1,[C(j),C(j-1)],sol1);
    [L2,R2]=AM3(dt,1,B2,[C(j+1),C(j),C(j-1)],sol2);
    [L3,R3]=BDF2(dt,1,B3,[C(j+1)],sol3);
    p4(j+1)=L1\R1;p5(j+1)=L2\R2;p6(j+1)=L3\R3;
    sol1=[p4(j+1),p4(j)];
    sol2=[p5(j+1),p5(j)];
    sol3=[p6(j+1),p6(j)];
end

plot(t,p1,'LineWidth',2),hold on
plot(t,p2,'LineWidth',2),hold on
plot(t,p3,'LineWidth',2),hold on
plot(t,p4,'LineWidth',2),hold on
plot(t,p5,'LineWidth',2),hold on
plot(t,p6,'LineWidth',2),hold on
plot(te,p,'k--','LineWidth',2),hold on
axis([0 2 0 5*10^-3])
legend('Euler Expl.','Euler Impl.','Trapez','AB2','AM3','BDF2','exakt')
 
 

 
常微分方程初值问题数值解法
 
问题
一、一阶常微分方程初值问题的有限差分方法与误差分析
二、向前Euler法及误差分析
1.向前Euler法
2.误差分析
3.后退Euler法
  
三、改进欧拉公式
四、单步法局部截断误差与阶
五、龙格—库塔方法
1.定义
2.常用的龙格库塔方法
  
六、单步法的收敛性与稳定性
1.收敛性与相容性
2.绝对稳定性与绝对稳定域
 

 
问题
 
一阶常微分方程的初值问题:
 $y^{{}^{\prime }}=f(x,y)$
 $y(x_0)=y_0$
 定理一
 设f在区域$D=(x,y)|a\le x\le b,y\in R$上连续，关于y满足$|f(x,y_1)-f(x,y_2)|\le L|y_1-y_2|,\forall x_0\in [a,b],y_0\in R$（这条件叫做Lipschitz条件）常微分方程的初值问题存在唯一的连续可微解$y(x)$
 
一、一阶常微分方程初值问题的有限差分方法与误差分析
 
只有一些特殊的常微分方程能求解析解，实际问题中主要靠数值解，所谓数值解法，就是寻求解y(x)在一系列离散节点$x_0<x_1<x_2......$上的近似值$y_0,y_1,y_2......$
 整体截断误差
 $g_i=|y_i-y(x_i)|$
 其中$y(x_i)$就是我们的精确解，$y_i$就是给定计算格式（迭代公式，初值是精确值）之后所计算出来的数值解。
 局部截断误差
 $e_i=|y(x_{i+1})-w_{i+1}|$
 其中$w_{i+1}$是 初值为$y(x_i)$ 时迭代出的数值解
 对于某种方法，局部截断误差满足$e_i=O(h^{p+1})$，则称该方法具有p阶精度。
 
二、向前Euler法及误差分析
 
1.向前Euler法
 
先设步长相等，大小为h,对常微分方程$y^{{}^{\prime }}=f(x,y)$中的倒数用均差近似，即
 $\frac{y(x_{i+1})-y(x_i)}{h}\approx y^{{}^{\prime }}=f(x_i,y(x_i))$
 若初值已知，利用迭代格式
 $$\{\begin{array}{l}{y}_0=y(x_0)\\ y_{i+1}=y_i+hf(x_i,y_i)\end{array}$$
 可逐次算出$y_1,y_2.....$
 
2.误差分析
 
对于精确值做泰勒展开式
 $y(x_{i+1})=y(x_i)+h{y}^{{}^{\prime }}(x_i)+\frac{h^2}{2}{y}^{{}^{\prime \prime }}(c)$其中$c\in (x_i,x_{i+1})$，其中$y^{{}^{\prime }}(x_i)=f(x_i,y(x_i))$
 对于欧拉法的迭代格式，由局部截断误差定义可得
 $y_{i+1}=y(x_i)+h{y}^{{}^{\prime }}(x_i)$
 局部截断误差$$e_i=\frac{h^2}{2}{y}^{{}^{\prime \prime }}(c)$$
 
3.后退Euler法
 
通过利用均差近似倒数$y^{{}^{\prime }}(x_{i+1})$即
 
$\frac{y(x_{i+1})-y(x_i)}{h}\approx y^{{}^{\prime }}(x_{i+1})=f(x_{i+1},y(x_{i+1}))$
 因此后退欧拉法的递推公式为
 $y_{i+1}=y_i+hf(x_{i+1},y_{i+1})$（1式）
 
相比欧拉法，后退欧拉法的递推公式并不能直接进行计算，而需要解方程，即此公式是隐式的，而欧拉法的递推公式是可以直接计算的，即为显式的。隐式方程通常用迭代法求解，而迭代过程的实质是逐步隐式化。
 $$\{\begin{array}{l}{y}_{i+1}^0=y_i+hf(x_i,y_i)\\ y_{i+1}^1=y_i+hf(x_i,y_{i+1}^0)\\ ......\\ y_{i+1}^{k+1}=y_i+hf(x_i,y_{i+1}^k)（2式）\end{array}$$
 由于f(x,y)满足Lipschitz条件，利用（1）式减去（2）式得到
 $|y_{i+1}^{k+1}-y_{i+1}|=h(f(x_i,y_{i+1}^k)-hf(x_{i+1},y_{i+1}))\le hL|y_{i+1}^k-y_{i+1}|$
 若$hL<1$则迭代法收敛到解$y_{i+1}$。
 
三、改进欧拉公式
 
补充：梯形方法
 $y_{i+1}=y_i+\frac{h}{2}(f(x_{i+1},y_{i+1})+f(x_i,y_i))$
 
先用欧拉法求得一个初步的近似值${\overline{y}}_{i+1}$，称为预报值，然后用它替代梯形法右端的$y_{i+1}$，得到左端校正值$y_{i+1}$,可以校正n次，但是改进欧拉法中，我们只校正一次，这样建立的预报-校正系统称为改进的欧拉格式:
 $$\{\begin{array}{l}{\overline{y}}_{i+1}=y_i+hf(x_i,y_i)\\ y_{i+1}=y_i+\frac{h}{2}(f(x_{i+1},{\overline{y}}_{i+1})+f(x_i,y_i))\end{array}$$
 也可以写成下面一个形式：
 $$\{\begin{array}{l}{y}_p=y_n+hf(x_n,y_n)\\ y_c=y_n+hf(x_n,y_p)\\ y_{n+1}=\frac{1}{2}(y_p+y_c)\end{array}$$
 
四、单步法局部截断误差与阶
 
初值问题$y^{{}^{\prime }}=f(x,y)$，$y(x_0)=y_0$的单步法可用一般形式表示为
 $y_{n+1}=y_n+h\phi (x_n，y_n，y_{n+1}，h)$
 其中$\phi (x_n，y_n，y_{n+1}，h)$称为增量函数
 显式单步法的表示形式
 $y_{n+1}=y_n+h\phi (x_n，y_n，h)$
 例如：对欧拉法$\phi (x_n，y_n，h)=f(x,y)$
 对于显式单步法的局部截断误差
 $T_{n+1}=y(x_{n+1})-y(x_n)-h\phi (x_n，y(x_n)，h)$
 
五、龙格—库塔方法
 
1.定义
 
对于方程$y^{{}^{\prime }}=f(x,y)$的等价积分形式
 
$y(x_{n+1})-y(x_n)={\int }_{x_n}^{x_{n+1}}f(x,y(x))dx$
 
方程右端可用求积公式表示为
 ${\int }_{x_n}^{x_{n+1}}f(x,y(x))dx\approx h{\sum }_{i=1}^r{c}_{i}f(x_n+{\lambda }_{i}h,y(x_n+{\lambda }_{i}h))$
 一般来说点数r越多精度越高，上式右端相当于增量函数$\phi (x_n，y_n，h)$，为了得到便于计算的显示方法，将公式表示为
 $y_{n+1}=y_n+h\phi (x_n,y_n,h)$
 其中
 $\phi (x_n,y_n,h)={\sum }_{i=1}^n{c}_i{K}_i$
 $K_1=f(x_n,y_n)$
 $K_2=f(x_n+{\lambda }_{2}h,y_n+h{u}_{21}{K}_1)，{\lambda }_2=h$
 …
 $K_i=f(x_n+{\lambda }_{i}h,y_n+h\times {\sum }_{j=1}^{i-1}{u}_{ij}{K}_j)，{\lambda }_i={\sum }_{j=1}^{i-1}{u}_{ij}$
 将此方法称为龙格库r级塔显示方法
 
2.常用的龙格库塔方法
 
常用二阶显式龙格库塔方法
 $y_{n+1}=y_n+h{K}_2$
 $K_1=f(x_n,y_n)$
 $K_2=f(x_n+\frac{h}{2},y_n+\frac{h}{2}{K}_1)$
 常用三阶
 $y_{n+1}=y_n+\frac{h}{6}(K_1+c2K_2)$
 $K_1=f(x_n,y_n)$
 $K_2=f(x_n+\frac{h}{2},y_n+\frac{h}{2}{K}_1)$
 $K_3=f(x_n+h,y_n-h{K}_1+2h{K}_2)$
 常用四阶（记下来，要考，哈哈哈,matlab计算内核所在）
 $y_{n+1}=y_n+\frac{h}{6}(K_1+2K_2+2K_3+K_4)$
 $K_1=f(x_n,y_n)$
 $K_2=f(x_n+\frac{h}{2},y_n+\frac{h}{2}{K}_1)$
 $K_3=f(x_n+\frac{h}{2},y_n+\frac{h}{2}{K}_2)$
 $K_4=f(x_n+h,y_n+h{K}_3)$
 
六、单步法的收敛性与稳定性
 
1.收敛性与相容性
 
定义：若一种数值方法（如单步法）对于固定的$x_n=x_0+nh$,当$h\to 0$ ，有$y_n=y(x_n)$,则称该方法收敛。
 定理：假设单步法具有p阶精度，且增量函数$\phi (x，y，h)$满足利普希茨条件
 $|\phi (x，y，h)-\phi (x，\overline{y}，h)|\le L_{\phi }|y-\overline{y}|$
 又设初值准确，则其整体截断误差
 $y(x_n)-y_n=O(h^p)$
 定义：若单步法$y_{n+1}=y_n+h\phi (x_n,y_n,h)$的增量函数满足$\phi (x,y,0)=f(x,y)$，则称单步法与初值问题相容（前述的单步法是用了泰勒展开式略去高阶项近似，有用求积公式离散得到的数值方法，相容性指数值方法逼近微分方程，在$h\to 0$时，可得到$y^{{}^{\prime }}(x)=f(x,y)$）
 
2.绝对稳定性与绝对稳定域
 
定义：若一种数值方法在节点值$y_n$上大小扰动为$\delta$，于以后的节点值$y_m(m>n)$上产生的偏差均不超过$\delta$，则称该方法是稳定的
 定义：对于模型方程$y^{{}^{\prime }}=\lambda y$,使用单步法解模型方程$y_{n+1}=E(\lambda h)y_n$，若得到$|E(h\lambda )|<1$，则称该方法是绝对稳定的（以此可以来确定步长）。
 常用：对于四阶龙格库塔方法:$-2.78<h\lambda <0$,即$0<h<-2.78/\lambda$

 
求解初值问题
 
简介
前期准备
欧拉法
改进的欧拉法
龙格-库塔法
标准四阶显式Kutta公式
三级三阶显式公式
四级四阶显式Kutta公式
四级四阶显式Gill公式
  
示例
MATLAB代码
结果
 

 
简介
 
通过求解简单的初值问题：
 $$\{\begin{array}{llccccc}\frac{du}{dx}=f(x,u)& & & & & & (1)\\ u(x_0)=u_0& & & & & & (2)\end{array}$$
 引入欧拉法、改进的欧拉法、龙格-库塔法等。
 
前期准备
 
数值解法的基本思想就是先对$x$和$u(x)$在区间$[x_0,\infty )$上进行离散化，然后构造递推公式，再进一步得到$u(x)$在这些位置的近似取值。
 
取定步长$h$，令$x_n=x_0+nh(n=\pm 1,\pm 2,\cdots )$
得到离散的位置：$x_1,x_2,\cdots ,x_n,\cdots$
$u(x)$在这些点精确取值为：$u(x_1),u(x_2),\cdots ,u(x_n),\cdots$
利用数值解法得到的这些点的近似取值，记为$u_1,u_2,\cdots ,u_n,\cdots$
 
欧拉法
 
欧拉法的核心就是将导数近似为差商。
 将导数近似为向前差商，则有：
 $${\frac{du}{dx}|}_{x=x_n}\approx \frac{u(x_{n+1})-u(x_n)}{h}$$
 代入(1)式，有：
 $$u(x_{n+1})=y(x_n)+hf(x_n,u(x_n))$$
 用$u_{n+1}$和$u_n$代替$u(x_{n+1})$和$u(x_n)$，得：
 $$u_{n+1}=u_n+hf(x_n,u_n)$$
 因此，若知道$u_0$我们就可以递归出$u_1,u_2,\cdots$
 如果将导数近似为向后差商:
 $${\frac{du}{dx}|}_{x=x_n}\approx \frac{u(x_n)-u(x_{n-1})}{h}$$
 类似的，就可以得到：
 $$u_{n-1}=u_n-hf(x_n,u_n)$$
 这样，若知道$u_0$我们就可以递归出$u_{-1},u_{-2},\cdots$
 
改进的欧拉法
 
对$(1)$式在$[x_n,x_{n+1}]$上积分，可得：
 $$u(x_{n+1})=u(x_n)+{\int }_{x_n}^{x_{n+1}}f(x,u)dx$$
 其中，$n=0,1,\cdots$。用不同方式来近似上式的积分运算，就会得到不同的递推公式。若使用左端点计算矩形面积并取近似：
 $${\int }_{x_n}^{x_{n+1}}f(x,u)dx\approx hf(x_{n+1},u(x_{n+1}))$$
 代入上式得：
 $$u_{n+1}=u_n+hf(x_n,u_n)$$
 若使用梯形的面积做近似：
 $${\int }_{x_n}^{x_{n+1}}f(x,y)dx\approx \frac{h}{2}[f(x_n,u(x_n))+f(x_{n+1},u(x_{n+1}))]$$
 得到：
 $$u_{n+1}=u_n+\frac{h}{2}[f(x_n,u_n)+f(x_{n+1},u_{n+1})]$$
 欧拉法虽然精度偏低，但它是显式的，可直接得到结果。而梯形公式是隐式的，虽然精度较高，却无法通过一步计算得到结果，若用迭代法计算，运算量较大。综合这两种方法，可以相得益彰：先用显式格式却低精度的欧拉法计算得到一个粗略的预测值${\bar{u}}_{n+1}$，再将这个预测值代入梯形公式进行修正，得到较高精度的结果$u_{n+1}$。
 $$\{\begin{array}{l}{\bar{u}}_{n+1}=u_n+hf(x_n,u_n)\\ u_{n+1}=u_n+\frac{h}{2}[f(x_n,u_n)+f(x_{n+1},{\bar{u}}_{n+1})]\end{array}$$
 
龙格-库塔法
 
将以上两种方法分别写成如下形式：
 $$\{\begin{array}{l}{u}_{n+1}=u_n+h{K}_1\\ K_1=f(x_n,u_n)\end{array}$$
 $$\{\begin{array}{l}{u}_{n+1}=u_n+\frac{h}{2}(K_1+K_2)\\ K_1=f(x_n,u_n)\\ K_2=f(x_n+h,u_n+K_1)\end{array}$$
 上述方法都是通过$f(x,u)$在不同位置的线性组合来计算$u_{n+1}$的值，所考虑的位置越多，精度也越高。类似的，就得到龙格-库塔法的思想：如果用$f(x,u)$在更多位置的线性组合来构造递推公式，将会得到更高的精度。
 这样，递推公式将有如下形式：
 $$\{\begin{array}{l}{u}_{n+1}=u_n+h\sum _{i=1}^r{R}_i{K}_i\\ K_1=f(x_n,u_n)\\ K_i=f(x_n+a_{i}h,u_n+\sum _{j=1}^{i-1}{b}_{ij}{K}_j),i=2,3,\cdots ,r\end{array}$$
 其中，$R_i,a_i,b_{ij}$为待定常数。(利用$Taylor$展开就可以确定待定系数)
 
标准四阶显式Kutta公式
 
$$\{\begin{array}{l}{y}_{n+1}=y_n+\frac{h}{6}(K_1+4K_2+K_3),\\ K_1=f(x_n,y_n),\\ K_2=f(x_n+\frac{1}{2}h,y_n+\frac{1}{2}h{K}_1),\\ K_3=f(x_n+h,y_n-h{K}_1+2h{K}_2);\end{array}$$
 
三级三阶显式公式
 
$$\{\begin{array}{l}{y}_{n+1}=y_n+\frac{h}{4}(K_1+3K_3),\\ K_1=f(x_n,y_n),\\ K_2=f(x_n+\frac{1}{3}h,y_n+\frac{1}{3}h{K}_1),\\ K_3=f(x_n+\frac{2}{3}h,y_n+\frac{2}{3}h{K}_2);\end{array}$$
 
四级四阶显式Kutta公式
 
$$\{\begin{array}{l}{y}_{n+1}=y_n+\frac{h}{8}(K_1+3K_2+3K_3+K_4),\\ K_1=f(x_n,y_n),\\ K_2=f(x_n+\frac{1}{3}h,y_n+\frac{1}{3}h{K}_1),\\ K_3=f(x_n+\frac{2}{3}h,y_n-\frac{2}{3}h{K}_1+h{K}_2),\\ K_4=f(x_n+h,y_n+h{K}_1-h{K}_2+h{K}_3);\end{array}$$
 
四级四阶显式Gill公式
 
$$\{\begin{array}{l}{y}_{n+1}=y_n+\frac{h}{6}(K_1+(2-\sqrt{2})K_2+(2+\sqrt{2})K_3+K_4),\\ K_1=f(x_n,y_n),\\ K_2=f(x_n+\frac{1}{2}h,y_n+\frac{1}{2}h{K}_1),\\ K_3=f(x_n+\frac{1}{2}h,y_n+\frac{\sqrt{2}-1}{2}h{K}_1+(1-\frac{\sqrt{2}}{2})h{K}_2),\\ K_4=f(x_n+h,y_n-\frac{\sqrt{2}}{2}h{K}_2+(1+\frac{\sqrt{2}}{2})h{K}_3);\end{array}$$
 
示例
 
求解常微分方程:
 $$\{\begin{array}{l}\frac{dy}{dx}=x^3-\frac{y}{x},\\ y(1)=\frac{2}{5}.\end{array}$$
 要求步长为$h=0.1$，其中，精确解为
 $$y=\frac{1}{5}{x}^4+\frac{1}{5x}.$$
 
MATLAB代码
 
clear all,close all,clc
f=@(x,y)x^3-y/x;
h=0.1;
%% Euler method
x=[1:h:2];
N=size(x,2)-1;
y1=[2/5,zeros(1,N)];
for n=1:N
    y1(n+1)=y1(n)+h*f(x(n),y1(n));
end
%% Improved Euler method
y2=[2/5,zeros(1,N)];
for n=1:N
    y2(n+1)=y2(n)+h*f(x(n),y2(n));
    y2(n+1)=y2(n)+h/2*(f(x(n),y2(n))+f(x(n+1),y2(n+1)));
end
%% Standard fourth-order explicit Kutta formula
y3=[2/5,zeros(1,N)];
for n=1:N
    K1=f(x(n),y3(n));
    K2=f(x(n)+1/2*h,y3(n)+1/2*h*K1);
    K3=f(x(n)+h,y3(n)-h*K1+2*h*K2);
    y3(n+1)=y3(n)+h/6*(K1+4*K2+K3);
end
%% Three-level three-order explicit formula
y4=[2/5,zeros(1,N)];
for n=1:N
    K1=f(x(n),y4(n));
    K2=f(x(n)+1/3*h,y4(n)+1/3*h*K1);
    K3=f(x(n)+2/3*h,y4(n)+2/3*h*K2);
    y4(n+1)=y4(n)+h/4*(K1+3*K3);
end
%% Fourth-level fourth-order explicit Kutta formula
y5=[2/5,zeros(1,N)];
for n=1:N
    K1=f(x(n),y5(n));
    K2=f(x(n)+1/3*h,y5(n)+1/3*h*K1);
    K3=f(x(n)+2/3*h,y5(n)-1/3*h*K1+h*K2);
    K4=f(x(n)+h,y5(n)+h*K1-h*K2+h*K3);
    y5(n+1)=y5(n)+h/8*(K1+3*K2+3*K3+K4);
end
%% Fourth-level fourth-order explicit Gill formula
y6=[2/5,zeros(1,N)];
for n=1:N
    K1=f(x(n),y6(n));
    K2=f(x(n)+1/2*h,y6(n)+1/2*h*K1);
    K3=f(x(n)+1/2*h,y6(n)+(sqrt(2)/2-0.5)*h*K1+(1-sqrt(2)/2)*h*K2);
    K4=f(x(n)+h,y6(n)-sqrt(2)/2*h*K2+(1+sqrt(2)/2)*h*K3);
    y6(n+1)=y6(n)+h/6*(K1+(2-sqrt(2))*K2+(2+sqrt(2))*K3+K4);
end
y=1/5*x.^4+1./(5*x);  %  Exact solution

plot(x,y,'k',x,y1,'xr',x,y2,'ob','Markersize',10,'LineWidth',1.5)
axis([1 2.1 0 4.5]),xlabel x,ylabel u
legend('Exact','Euler','Improved Euler')

%plot(x,y,'k',x,y3,'*r',x,y4,'+b',x,y5,'-y',x,y6,'or','Markersize',10,'LineWidth',1.5)
%axis([1 2.1 0 6]),xlabel x,ylabel u,set(gca,'Fontsize',18)
%legend('Exact','34Kutta','33kutta','44Kutta','44Gill')

 
结果
 

介绍
 
关于初值问题的解
 
 
 
实例
 
 
代码
 
#include<iostream>
using namespace std;

double h;
double x[100],y[100];

double f(double x0, double y0)
{
	return -2*x0/y0 + y0;	
}

void Eluer()
{
	for(int i=1;i<10;i++)
	{
		y[i] = y[i-1] + h*f(x[i-1],y[i-1]);
		cout<<x[i]<<" "<<y[i]<<endl;
	}
}

int main()
{
	cout<<"输入迭代的步长："<<endl;
	cin>>h;
	cout<<"输入初始值x,y"<<endl;
	cin>>x[0]>>y[0];
	for(int i=1;i<10;i++)  x[i] = x[0] + i*h;
//	for(int i=1;i<10;i++)  cout<<x[i]<<" ";cout<<endl; 
	Eluer();
	return 0;
}
 
结果(解与实例中对应 ):
 
 
Euler 改进算法
 
 
实例
 
 
代码
 
#include<iostream>
using namespace std;

double h;
double x[100],y[100],c[100];

double f(double x0, double y0)
{
	return -2*x0/y0 + y0;	
}

void Eluer()
{
	for(int i=1;i<10;i++)
	{
		c[i] = y[i-1] + h*f(x[i-1],y[i-1]);
		y[i] = y[i-1] + h/2*(f(x[i-1],y[i-1])+f(x[i],c[i]));
		cout<<x[i]<<" "<<y[i]<<endl;
	}
}

int main()
{
	cout<<"输入迭代的步长："<<endl;
	cin>>h;
	cout<<"输入初始值x,y"<<endl;
	cin>>x[0]>>y[0];
	for(int i=1;i<10;i++)  x[i] = x[0] + i*h;
//	for(int i=1;i<10;i++)  cout<<x[i]<<" ";cout<<endl; 
	Eluer();
	return 0;
}
 
结果(解与实例中对应 ):