Matlab实验报告五(微分方程求解Euler折线法)-推荐下载

实验五 利用Matlab 求解常微分方程(组)的实验报告学院：数计学院 班级：1003班 姓名：黄晓丹 学号：1051020144一.实验目的:熟悉Matlab 软件中关于求解常微分方程的各种命令. 掌握利用Matlab 软件进行常微分方程的求解。

二.相关知识在MATLAB 中，由函数dsolve()解决常微分方程（组）的求解问题，其具体格式如下：X=dsolve(‘eqn1’,’eqn2’,…)函数dsolve 用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解.用字符串方程表示，自变量缺省值为t 。

导数用D 表示，2阶导数用D2表示，以此类推。

S 返回解析解。

在方程组情形，s 为一个符号结构。

[tout,yout]=ode45(‘yprime’,[t0,tf],y0) 采用变步长四阶Runge-Kutta 法和五阶Runge-Kutta-Felhberg 法求数值解，yprime 是用以表示f(t,y)的M 文件名，t0表示自变量的初始值，tf 表示自变量的终值，y0表示初始向量值。

输出向量tout 表示节点(t 0,t 1, …,t n )T ,输出矩阵yout 表示数值解，每一列对应y 的一个分量。

若无输出参数，则自动作出图形.三.实验内容:例1 求下列微分方程的解析解（1）（2）（3）方程（1）求解的MATLAB 代码为：>>clear;>>s=dsolve('Dy=a*y+b')结果为s =-b/a+exp(a*t)*C1方程（2）求解的MATLAB 代码为：>>clear;>>s=dsolve('D2y=sin(2*x)-y','y(0)=0','Dy(0)=1','x')>>simplify(s) %以最简形式显示s结果为s =(-1/6*cos(3*x)-1/2*cos(x))*sin(x)+(-1/2*sin(x)+1/6*sin(3*x))*cos(x)+5/3*sin(x)ans =-2/3*sin(x)*cos(x)+5/3*sin(x)方程（3）求解的MATLAB 代码为：b ay y +='1)0(',0)0(,)2sin(''==-=y y y x y 1)0(',1)0(',','==-=+=g f f g g g f f>>clear;>>s=dsolve('Df=f+g','Dg=g-f','f(0)=1','g(0)=1')>>simplify(s,f) %s 是一个结构>>simplify(s.g)结果为ans =exp(t)*cos(t)+exp(t)*sin(t)ans =-exp(t)*sin(t)+exp(t)*cos(t)例２ 求解微分方程先求解析解，再求数值解，并进行比较。

matlab欧拉法求解微分方程例题当使用Matlab的欧拉法（Euler's method）求解微分方程时，需要将微分方程转化为离散的差分方程。

下面以一个简单的一阶常微分方程为例来说明。

假设我们要求解以下的微分方程：dy/dx = x + y并给定初始条件 y(0) = 1。

我们可以使用欧拉法来逼近该微分方程的解。

在Matlab中，可以按照以下步骤进行求解：1. 定义步长和计算步数：h = 0.1; % 步长 N = 10; % 计算步数2. 初始化变量并给定初始条件：x = zeros(N+1, 1); % 初始化 x 数组 y = zeros(N+1, 1); % 初始化 y 数组 x(1) = 0; % 初始条件 x(0) = 0 y(1) = 1; % 初始条件 y(0) = 13. 使用欧拉法进行迭代计算：for i = 1:N x(i+1) = x(i) + h; % 计算下一个 x 值 y(i+1) = y(i) + h * (x(i) + y(i)); % 根据差分方程计算下一个 y 值end4. 绘制结果：plot(x, y, 'o-'); % 绘制连线图 xlabel('x'); ylabel('y'); title('Approximate Solution using Euler''s Method');完整的代码如下：h = 0.1; % 步长 N = 10; % 计算步数 x = zeros(N+1, 1); % 初始化 x 数组 y = zeros(N+1, 1); % 初始化 y 数组 x(1) = 0; % 初始条件 x(0) = 0 y(1) = 1; % 初始条件 y(0) = 1 for i = 1:N x(i+1) = x(i) + h; % 计算下一个 x 值y(i+1) = y(i) + h * (x(i) + y(i)); % 根据差分方程计算下一个 y 值 end plot(x, y, 'o-'); % 绘制连线图xlabel('x'); ylabel('y'); title('Approximate Solution using Euler''s Method');运行该代码，将得到近似的解曲线图，表示微分方程的数值解。

实验二微分方程求解一、问题背景与实验目的实际应用问题通过数学建模所归纳而取得的方程，绝大多数都是微分方程，真正能取得代数方程的机遇很少．另一方面，能够求解的微分方程也是十分有限的，专门是高阶方程和偏微分方程（组）．这就要求咱们必需研究微分方程（组）的解法，既要研究微分方程（组）的解析解法（精准解），更要研究微分方程（组）的数值解法（近似解）．对微分方程（组）的解析解法(精准解)，Matlab 有专门的函数能够用，本实验将作必然的介绍．本实验将要紧研究微分方程(组)的数值解法（近似解），重点介绍 Euler 折线法．二、相关函数（命令）及简介1．dsolve('equ1','equ2',…)：Matlab 求微分方程的解析解．equ一、equ 二、…为方程（或条件）．写方程（或条件）时用 Dy 表示y 关于自变量的一阶导数，用用 D2y 表示 y 关于自变量的二阶导数，依此类推．2．simplify(s)：对表达式 s 利用 maple 的化简规那么进行化简．例如：syms xsimplify(sin(x)^2 + cos(x)^2)ans=13．[r,how]=simple(s)：由于 Matlab 提供了多种化简规那么，simple 命令确实是对表达式 s 用各类规那么进行化简，然后用 r 返回最简形式，how 返回形成这种形式所用的规那么．例如：syms x[r,how]=simple(cos(x)^2-sin(x)^2)r = cos(2*x)how = combine4．[T,Y] = solver(odefun,tspan,y 0) 求微分方程的数值解．说明：(1) 其中的 solver 为命令 ode4五、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 之一．(2) odefun 是显式常微分方程：⎪⎩⎪⎨⎧==00)(),(y t y y t f dt dy(3) 在积分区间 tspan =],[0f t t 上，从0t 到f t ，用初始条件0y 求解．(4) 要取得问题在其他指按时刻点 ,210,,t t t 上的解，那么令 tspan =],,,[,210f t t t t （要求是单调的）． (5) 因为没有一种算法能够有效地解决所有的 ODE 问题，为此，Matlab 提供了多种求解器 Solver ，关于不同的ODE 问题，采纳不同的Solver ．(6) 要专门的是：ode23、ode45 是极为经常使用的用来求解非刚性的标准形式的一阶常微分方程(组)的初值问题的解的 Matlab 的经常使用程序，其中：ode23 采纳龙格-库塔2 阶算法，用3 阶公式作误差估量来调剂步长，具有低等的精度．ode45 那么采纳龙格-库塔4 阶算法，用5 阶公式作误差估量来调剂步长，具有中等的精度．5．ezplot(x,y,[tmin,tmax])：符号函数的作图命令．x,y 为关于参数t 的符号函数，[tmin,tmax] 为 t 的取值范围．6．inline()：成立一个内联函数．格式：inline('expr', 'var1', 'var2',…) ，注意括号里的表达式要加引号．例：Q = dblquad(inline('y*sin(x)'), pi, 2*pi, 0, pi)三、实验内容1. 几个能够直接用 Matlab 求微分方程精准解的例子：例1：求解微分方程22x xe xy dxdy -=+，并加以验证． 求解本问题的Matlab 程序为：syms x y %line1y=dsolve('Dy+2*x*y=x*exp(-x^2)','x') %line2diff(y,x)+2*x*y-x*exp(-x^2) %line3simplify(diff(y,x)+2*x*y-x*exp(-x^2)) %line4说明：(1) 行line1是用命令概念x,y 为符号变量．那个地址能够不写，但为确保正确性，建议写上；(2) 行line2是用命令求出的微分方程的解：1/2*exp(-x^2)*x^2+exp(-x^2)*C1(3) 行line3利用所求得的解．那个地址是将解代入原微分方程，结果应该为0，但那个地址给出：-x^3*exp(-x^2)-2*x*exp(-x^2)*C1+2*x*(1/2*exp(-x^2)*x^2+exp(-x^2)*C1)(4) 行line4 用 simplify() 函数对上式进行化简，结果为 0， 说明)(x y y =的确是微分方程的解．例2：求微分方程0'=-+x e y xy 在初始条件e y 2)1(=下的特解，并画出解函数的图形．求解本问题的 Matlab 程序为：syms x yy=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x')ezplot(y)微分方程的特解为：y=1/x*exp(x)+1/x* exp (1) (Matlab 格式)，即xe e y x+=，解函数的图形如图 1：图1例3：求微分方程组⎪⎪⎩⎪⎪⎨⎧=--=++035y x dtdy e y x dt dx t 在初始条件0|,1|00====t t y x 下的特解，并画出解函数的图形．求解本问题的 Matlab 程序为：syms x y t[x,y]=dsolve('Dx+5*x+y=exp(t)','Dy-x-3*y=0','x(0)=1','y(0)=0','t')simple(x);simple(y);ezplot(x,y,[0,]);axis auto微分方程的特解(式子专门长)和解函数的图形均略．2. 用ode23、ode45等求解非刚性的标准形式的一阶常微分方程(组)的初值问题的数值解(近似解)．例4：求解微分方程初值问题⎪⎩⎪⎨⎧=++-=1)0(2222y x x y dx dy 的数值解，求解范围为区间[0, ]．fun=inline('-2*y+2*x^2+2*x','x','y');[x,y]=ode23(fun,[0,],1);x';y';plot(x,y,'o-')>> x'ans =>> y'ans =图形结果为图 2．图2例 5：求解描述振荡器的经典的 Ver der Pol 微分方程.7,0)0(',1)0(,0)1(222====+--μμy y y dt dy y dty d 分析：令,,121dt dx x y x ==则.)1(,1221221x x x dtdx x dt dx --==μ 先编写函数文件： function xprime = verderpol(t,x)global mu;xprime = [x(2);mu*(1-x(1)^2)*x(2)-x(1)];再编写命令文件：global mu;mu = 7;y0=[1;0][t,x] = ode45('verderpol',[0,40],y0);x1=x(:,1);x2=x(:,2);plot(t,x1)图形结果为图3．图33. 用 Euler 折线法求解前面讲到过，能够求解的微分方程也是十分有限的．下面介绍用 Euler 折线法求微分方程的数值解(近似解)的方式．Euler 折线法求解的大体思想是将微分方程初值问题⎪⎩⎪⎨⎧==00)(),,(y x y y x f dx dy 化成一个代数方程，即差分方程，要紧步骤是用差商h x y h x y )()(-+替代微商dxdy ,于是： ⎪⎩⎪⎨⎧==-+)()),(,()()(00x y y x y x f h x y h x y k k k k 记)(,1k k k k x y y h x x =+=+，从而)(1h x y y k k +=+，那么有1,,2,1,0).,(,),(1100-=⎪⎩⎪⎨⎧+=+==++n k y x hf y y h x x x y y k k k k k k 例 6：用 Euler 折线法求解微分方程初值问题⎪⎩⎪⎨⎧=+=1)0(,22y y x y dx dy的数值解(步长h 取，求解范围为区间[0,2]．解：本问题的差分方程为1,,2,1,0).2),( ),(,,4.0,1,021100-=⎪⎪⎪⎩⎪⎪⎪⎨⎧+=+=+====++n k y x y y x f y x hf y y h x x h y x k k k k k k （其中： 相应的Matlab 程序见附录 1．附录 1：clearf=sym ('y+2*x/y^2');a=0;b=2;h=;n=(b-a)/h+1;x=0;y=1;szj=[x,y];for i=1:n-1y=y+h*subs(f,{'x','y'},{x,y});x=x+h;szj=[szj;x,y];endszjplot(szj(:,1),szj(:,2))数据结果为：图形结果见图4：图4专门说明：本问题可进一步利用四阶 Runge-Kutta 法求解，读者可将两个结果在一个图中显示，并和精准值比较，看看哪个更“精准”？（相应的 Matlab 程序参见附录 2）．四、自己动手1. 求微分方程0-x-yx的通解．xy+sin2')1(2=2. 求微分方程x2''='+-的通解．ye5yy x sin3. 求微分方程组⎪⎪⎩⎪⎪⎨⎧=-+=++00y x dtdy y x dt dx 在初始条件0|,1|00====t t y x 下的特解，并画出解函数()y f x =的图形．4. 别离用 ode23、ode45 求上述第 3 题中的微分方程初值问题的数值解(近似解)，求解区间为[0,2]t ∈．利用画图来比较两种求解器之间的不同．5. 用 Euler 折线法求解微分方程初值问题⎪⎩⎪⎨⎧=-=1)0(,12'32y y x y y 的数值解(步长h 取，求解范围为区间[0,2]．6. 用四阶 Runge-Kutta 法求解微分方程初值问题⎩⎨⎧=-=1)0(,cos 'y x e y y x 的数值解(步长h 取，求解范围为区间[0,3]．四阶 Runge-Kutta 法的迭代公式为(Euler 折线法实为一阶 Runge-Kutta 法)：1,,2,1,0),()2,2()2,2(),()22(6,),(342312143211100-=⎪⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎪⎨⎧++=++=++==++++=+==++n k hL y h x f L L h y h x f L L h y h x f L y x f L L L L L h y y h x x x y y k k k k k k k k k k k k 相应的 Matlab 程序参见附录 2．附录 2：clearf=sym('y-exp(x)*cos(x)');a=0;b=3;h=;n=(b-a)/h+1;x=0;y=1;szj=[x,y];for i=1:n-1l1=subs(f,{'x','y'},{x,y});l2=subs(f,{'x','y'},{x+h/2,y+l1*h/2});l3=subs(f,{'x','y'},{x+h/2,y+l2*h/2});l4=subs(f,{'x','y'},{x+h,y+l3*h});y=y+h*(l1+2*l2+2*l3+l4)/6;x=x+h;szj=[szj;x,y];endszjplot(szj(:,1),szj(:,2))试用该方式求解第5题中的初值问题．7. 用 ode45 方式求上述第 6 题的常微分方程初值问题的数值解(近似解)，从而利用画图来比较二者间的不同．五、附录附录 1：clearf=sym('y+2*x/y^2');a=0;b=2;h=;n=(b-a)/h+1;x=0;y=1;szj=[x,y];for i=1:n-1y=y+h*subs(f,{'x','y'},{x,y});x=x+h;szj=[szj;x,y];endszjplot(szj(:,1),szj(:,2))附录 2：clearf=sym('y-exp(x)*cos(x)');a=0;b=3;h=;n=(b-a)/h+1;x=0;y=1;szj=[x,y];for i=1:n-1l1=subs(f,{'x','y'},{x,y});l2=subs(f,{'x','y'},{x+h/2,y+l1*h/2});l3=subs(f,{'x','y'},{x+h/2,y+l2*h/2});l4=subs(f,{'x','y'},{x+h,y+l3*h});y=y+h*(l1+2*l2+2*l3+l4)/6;x=x+h;szj=[szj;x,y];endszjplot(szj(:,1),szj(:,2))下面是个例子：% 欲解方程组：% dx/dt=y-x+x(t-1);% dy/dt=y-xy+y(t-1);x0=-1; %初试条件y0=1; %初试条件t=0:01:3;yy=zeros(size(t));xx=zeros(size(t));yy(1)=x0;yy(1)=y0;dt=t(2)-t(1);Nk=round(1/dt);for k=1:length(t)-1;if t(k)<1;k1=yy(k)-xx(k);g1=yy(k)-xx(k)*yy(k);k2=yy(k)+*g1*dt-[xx(k)+*dt*k1];g2=yy(k)+*g1*dt-[xx(k)+*dt*k1]*[yy(k)+*g1*dt];k3=yy(k)+*g2*dt-[xx(k)+*dt*k2];g3=yy(k)+*g2*dt-[xx(k)+*dt*k1]*[yy(k)+*g2*dt];k4=yy(k)+g3*dt-[xx(k)+dt*k3];g4=yy(k)+g3*dt-[xx(k)+dt*k3]*[yy(k)+g3*dt];elsek1=yy(k)-xx(k)+xx(k-Nk);g1=yy(k)-xx(k)*yy(k)+yy(k-Nk);k2=yy(k)+*g1*dt-[xx(k)+*dt*k1]+xx(k-Nk)+*dt*k1;g2=yy(k)+*g1*dt-[xx(k)+*dt*k1]*[yy(k)+*g1*dt]+yy(k-Nk)+*dt*g1;k3=yy(k)+*g2*dt-[xx(k)+*dt*k2]+xx(k-Nk)+*dt*k2;g3=yy(k)+*g2*dt-[xx(k)+*dt*k1]*[yy(k)+*g2*dt]+yy(k-Nk)+*dt*g2;k4=yy(k)+g3*dt-[xx(k)+dt*k3]+xx(k-Nk)+dt*k3;g4=yy(k)+g3*dt-[xx(k)+dt*k3]*[yy(k)+g3*dt]+yy(k-Nk)+dt*k3;endxx(k+1)=xx(k)+dt*[k1+2*k2+2*k3+k4]/6;yy(k+1)=yy(k)+dt*[g1+2*g2+2*g3+g4]/6;endsubplot(121);plot(t,xx);title('t~x(t)');subplot(122);plot(t,yy);title('t~y(t)');/ / / / 文件夹:一、% ４阶龙格库塔法求chen's吸引子% 方程表达式% dx/dt = a*(y-x)% dy/dt = (c-a)*x - x*z + c*y% dz/dt = x*y - b*zclc;close all;clear all;%参数值a = 35;b = 3;c = 28;%初始值x_0 = -1;y_0 = 0;z_0 = 1;h = ; % 积分时间步长step1 = 10000; % 前面的迭代点数step2 = 5000; % 后面的迭代点数X = [];Y = [];Z = [];for(i = 1:1:(step1 + step2))%e1x_e1 = -a*x_0 + a*y_0;y_e1 = c*y_0 + (c - a)*x_0 - x_0*z_0;z_e1 = -b*z_0 + x_0*y_0;%e2y_h = y_0 + *h*y_e1;x_e2 = -a*(x_0 + *h*x_e1) + a*y_h;x_h = x_0 + *h*x_e1;z_h = z_0 + *h*z_e1;y_e2 = c*(y_0 + *h*y_e1) + (c - a)*x_h - x_h*z_h;z_e2 = -b*(z_0 + *h*z_e1) + x_h*y_h;%e3y_h = y_0 + *h*y_e2;x_e3 = -a*(x_0 + *h*x_e2) + a*y_h;x_h = x_0 + *h*x_e2;z_h = z_0 + *h*z_e2;y_e3 = c*(y_0 + *h*y_e2) + (c - a)*x_h - x_h*z_h;z_e3 = -b*(z_0 + *h*z_e2) + x_h*y_h;%e4y_h = y_0 + h*y_e3;x_e4 = -a*(x_0 + h*x_e3) + a*y_h;x_h = x_0 + h*x_e3;z_h = z_0 + h*z_e3;y_e4 = c*(y_0 + h*y_e2) + (c - a)*x_h - x_h*z_h;z_e4 = -b*(z_0 + h*z_e2) + x_h*y_h;%叠代x_1 = x_0 + 1/6*h*(x_e1 + 2*x_e2 +2*x_e3 + x_e4);y_1 = y_0 + 1/6*h*(y_e1 + 2*y_e2 +2*y_e3 + y_e4);z_1 = z_0 + 1/6*h*(z_e1 + 2*z_e2 +2*z_e3 + z_e4);X = [X,x_1];Y = [Y,y_1];Z = [Z,z_1];x_0 = x_1;y_0 = y_1;z_0 = z_1;endX = X(step2+1:end);Y = Y(step2+1:end);Z = Z(step2+1:end);figure(1);plot3(Z(1000:end),Y(1000:end),X(1000:end));grid onxlabel('x');ylabel('y');zlabel('z');这个应该是可以借鉴的。

MATLAB方法实现信号波形和微分方程的建立与求解实验报告实验目的：1.熟悉MATLAB软件的基本使用方法；2.掌握信号波形的绘制；3.掌握微分方程的建立与求解。

实验仪器：-计算机-MATLAB软件实验原理：1.信号波形的绘制信号波形是指随时间变化的量的图示。

可以使用MATLAB软件实现信号波形的绘制。

MATLAB软件提供了丰富的绘图函数，如plot函数、stem 函数等，可以根据需要选择合适的函数进行波形绘制。

2.微分方程的建立与求解微分方程是描述自然现象中变化的数学模型。

可以通过MATLAB软件进行微分方程的建立与求解。

MATLAB提供了ode45函数等函数用于求解常微分方程。

实验步骤：1.信号波形的绘制首先，打开MATLAB软件，并创建一个新的脚本文件。

然后，在脚本文件中使用plot函数绘制一个正弦信号的波形。

代码示例如下：```t=0:0.1:10;%时间范围为0到10，步长为0.1y = sin(t); % 正弦信号plot(t, y); % 绘制波形图```保存脚本文件，点击“运行”按钮即可显示出正弦信号的波形图。

2.微分方程的建立与求解在同一个脚本文件中，可以进行微分方程的建立与求解。

首先，使用syms命令定义微分方程中的未知函数和其导数。

然后，使用dsolve命令求解微分方程。

最后，通过plot函数绘制微分方程的解的波形。

代码示例如下：```syms y(t)ode = diff(y, t) == y; % 定义微分方程ySol(t) = dsolve(ode); % 求解微分方程t=0:0.1:10;%时间范围为0到10，步长为0.1y = ySol(t); % 微分方程的解plot(t, y); % 绘制波形图```保存脚本文件，点击“运行”按钮即可显示出微分方程的解的波形图。

实验结果：根据实验步骤，在MATLAB软件中实现了信号波形的绘制和微分方程的建立与求解。

通过绘图函数，成功绘制出了正弦信号的波形和微分方程的解的波形。

欧拉法是数值分析中常用的一种方法，用于求解常微分方程的数值解。

在MATLAB中，可以通过编写相应的代码来实现欧拉法求解微分方程。

下面我们将通过具体的实例来讲解MATLAB中如何使用欧拉法求解微分方程。

我们要了解欧拉法的基本原理。

欧拉法是一种通过迭代逼近微分方程解的方法，它基于微分方程的定义，通过离散化的方法逼近微分方程的解。

其基本思想是利用微分方程的导数定义，将微分方程以差分形式进行逼近。

具体而言，欧拉法通过将微分方程转化为差分方程的形式，然后通过迭代逼近得到微分方程的数值解。

接下来，我们通过一个具体的实例来讲解MATLAB中如何使用欧拉法求解微分方程。

假设我们要求解以下的一阶常微分方程：(1) dy/dx = x + y(2) y(0) = 1现在我们来编写MATLAB代码来实现欧拉法求解这个微分方程。

我们需要确定微分方程的迭代步长和迭代范围。

假设我们将x的范围取为0到10，步长为0.1。

接下来，我们可以编写MATLAB代码如下：```matlab欧拉法求解微分方程 dy/dx = x + y定义迭代步长和范围h = 0.1;x = 0:h:10;初始化y值y = zeros(1,length(x));y(1) = 1;使用欧拉法迭代求解for i = 1:(length(x)-1)y(i+1) = y(i) + h * (x(i) + y(i));end绘制图像plot(x,y,'-o');xlabel('x');ylabel('y');title('欧拉法求解微分方程 dy/dx = x + y');```在这段MATLAB代码中，我们首先定义了迭代的步长和范围，并初始化了微分方程的初始值y(0) = 1。

然后通过for循环使用欧拉法进行迭代求解微分方程，最后绘制出了微分方程的数值解的图像。

通过以上的实例讲解，我们可以看到，在MATLAB中使用欧拉法求解微分方程是非常简单而直观的。

佛山科学技术学院实 验 报 告课程名称 数值分析 实验项目 常微分方程问题初值问题数值解法 专业班级 姓 名 学 号 指导教师 陈剑 成 绩 日 期一. 实验目的1、理解如何在计算机上实现用Euler 法、改进Euler 法、Runge －Kutta 算法求一阶常微分方程初值问题⎩⎨⎧=∈='1)(],[),,()(y a y b a x y x f x y 的数值解。

2、利用图形直观分析近似解和准确解之间的误差。

二、实验要求（1） 按照题目要求完成实验内容； （2） 写出相应的Matlab 程序；（3） 给出实验结果(可以用表格展示实验结果)； （4） 分析和讨论实验结果并提出可能的优化实验。

（5） 写出实验报告。

三、实验步骤1、用Matlab 编写解常微分方程初值问题的Euler 法、改进Euler 法和经典的四阶Runge-Kutta 法。

2、给定初值问题⎪⎩⎪⎨⎧=≤≤-=;1)1(,21,1')1(2y x xy x y⎪⎩⎪⎨⎧=≤≤++-=31)0(10,25050')2(2y x x x y y 要求：（a ）用Euler 法和改进的Euler 法(步长均取h=0.05)及经典的四阶Runge-Kutta 法(h=0.1)求（1）的数值解，并打印)10,....2,1,0(1.01=+=i i x 的值。

(b) 用经典的四阶Runge-Kutta 方法解(2)，步长分别取h=0.1, 0.05,0.025,计算并打印)10,....2,1,0(1.0==i i x 个点的值，与准确解25031)(x e x y x +=-比较，并列表写出在x=0.2,0.5,0.8处，对于不同步长h 下的误差，讨论同一节点处，误差随步长的变化规律。

（c ）用Matlab 绘图函数绘制（2）的精确解和近似解的图形。

四、实验结果 %Euler.mfunction y = Euler(x0,xn,y0,h) %Euler 法解方程f_xy ； %x0,y0为初始条件； %x0,xn 为求值区间； %h 为步长； %求区间个数： n = (xn-x0)/h;%矩阵x 存储n+1个节点： x = [x0:h:xn]';%矩阵y 存储节点处的值： y = [y0;zeros(n,1)];%矩阵y_存储节点处导数值： y_(1)= f_xy(x0,y0); y_ = [y_(1);zeros(n,1)];%进行迭代（欧拉法迭代；求导数）： for i = 2:n+1y (i) = y(i-1)+h*y_(i-1); y_(i) = f_xy(x(i),y(i)); end%Imp_Euler.mfunction y = Imp_Euler(x0,xn,y0,h)%改进的Euler法解方程f_xy；%x0,y0为初始条件；%x0,xn为求值区间；%h为步长；%求区间个数：n = (xn-x0)/h;%矩阵x存储n+1个节点：x = [x0:h:xn]';%矩阵y存储节点处的值：y = [y0;zeros(n,1)];%矩阵y_存储节点处导数值：y_(1)= f_xy(x0,y0);y_ = [y_(1);zeros(n,1)];%使用改进Euler法求值（欧拉法求近似；近似点导数；梯形校正；求导）：for i = 2:n+1y_l = y(i-1) + h*y_(i-1);y_l = f_xy(x(i),y_l);y(i) = y(i-1) + (h/2)*(y_(i-1)+y_l);y_(i)= f_xy(x(i),y(i));end%R_Kutta4.mfunction y = R_Kutta4(x0,xn,y0,h)%Runger_Kutta法解方程f_xy；%x0,y0为初始条件；%x0,xn为求值区间；%h为步长；%求区间个数：n = (xn-x0)/h;%矩阵x存储n+1个节点：x = [x0:h:xn]';%矩阵y存储节点处的值：y = [y0;zeros(n,1)];%矩阵k1,k2,k3,k4存储各节点（中点）数值：k1(1)= f_xy(x0,y0);k1 = [k1(1);zeros(n,1)];k2(1)= f_xy(x0+h/2,y0+h*k1(1)/2);k2 = [k2(1);zeros(n,1)];k3(1)= f_xy(x0+h/2,y0+h*k2(1)/2);k3 = [k3(1);zeros(n,1)];k4(1)= f_xy(x0+h,y0+h*k3(1));k4 = [k4(1);zeros(n,1)];for i= 2:n+1y(i) = y(i-1)+(h/6)*(k1(i-1)+2*k2(i-1)+2*k3(i-1)+k4(i-1));k1(i)= f_xy(x(i),y(i));k2(i)= f_xy(x(i)+h/2,y(i)+h*k1(i)/2);k3(i)= f_xy(x(i)+h/2,y(i)+h*k2(i)/2);k4(i)= f_xy(x(i)+h,y(i)+h*k3(i));end（a）：%f_xy.mfunction y_=f_xy(x,y)%求解第五次作业第一题的点（x,y）处的导数；y_ = 1/(x^2) - y/x;%run521.mclc,clear;x0 = 1;xn = 2;h = 0.05;y0 = 1;%便于显示出x，与y对应：x = [x0:h:xn]';y = Euler(x0,xn,y0,h);YE =[x,y];y = Imp_Euler(x0,xn,y0,h); YIE= [x,y];h = 0.1;x = [x0:h:xn]';y = R_Kutta4(x0,xn,y0,h); YRK= [x,y];(b)： %f_xy.mfunction y_=f_xy(x,y) %求第二个方程的导数： y_ = -50*y+50*(x^2)+2*x;%run522.mclc,clear; x0 = 0; xn = 1; y0 = 1/3; %步长0.1： h = 0.1; x = [x0:h:xn]';y = R_Kutta4(x0,xn,y0,h); y_r= f_Real(x); Y1 = [x,y,y_r];%步长0.05： h = 0.05; x = [x0:h:xn]';y = R_Kutta4(x0,xn,y0,h); y_r= f_Real(x); Y2 = [x,y,y_r]; %步长0.025： h = 0.025; x = [x0:h:xn]';y = R_Kutta4(x0,xn,y0,h); y_r= f_Real(x); Y3 = [x,y,y_r];五、讨论分析(a)从结果可以看出使用RK 方法，步长较大但是结果也更加精确； (b)分析求值结果的误差，可以发现当步长取0.1时，误差是超级大的（10^8数量级），但是当步长缩小一半取0.05时，误差就很小了，再缩小一半，误差就更小了。