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

合集下载

用MATLAB求解常微分方程实验报告

用MATLAB求解常微分方程实验报告

实验五 利用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)例2 求解微分方程先求解析解,再求数值解,并进行比较。

matlab欧拉法求解微分方程例题

matlab欧拉法求解微分方程例题

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求微分方程的解

MATLAB求微分方程的解

实验二微分方程求解一、问题背景与实验目的实际应用问题通过数学建模所归纳而取得的方程,绝大多数都是微分方程,真正能取得代数方程的机遇很少.另一方面,能够求解的微分方程也是十分有限的,专门是高阶方程和偏微分方程(组).这就要求咱们必需研究微分方程(组)的解法,既要研究微分方程(组)的解析解法(精准解),更要研究微分方程(组)的数值解法(近似解).对微分方程(组)的解析解法(精准解),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)');/ / / / 文件夹:一、% 4阶龙格库塔法求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方法实现信号波形和微分方程的建立与求解实验报告

MATLAB方法实现信号波形和微分方程的建立与求解实验报告

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求微分方程精确解与近似解

Euler 折线法
考虑一维经典初值问题
dy dx
f (x, y) ,
y( x0 )
y0
,
x [a, b]
基本思想:用差商代替微商
根据 Talyor 公式,y(x) 在点 xk 处有
y( x) y( xk ) ( x xk ) y '( xk ) O(x2 )
y( xk1 ) y( xk ) hy '( xk ) O(h2 ) h xk1 xk
clear; global mu;mu=7; y0=[1;0]; [t,x]=ode45('verderpol',[0,40],y0); plot(t,x(:,1),'r-');
Matlab 求解微分方程小结
Matlab 函数
求解析解(通解或特解),用 dsolve 求数值解(特解),用 ode45、ode23 ...
Matlab 编程
Euler 折线法 Runga-Kutta 方法
上机作业
教材 P97:练习 1、2、3、4、5、6、7 要求:
请在上机之前将程序写好(参考附录),上机时 直接输入或修改附录程序即可。

4:求初值问题
dy dx
2 y 2x2
2x
的数值解,求解范
围为 [0,0.5]
y(0) 1
fun=inline('-2*y+2*x^2+2*x','x','y'); [x,y]=ode23(fun,[0,0.5],1);
注:也可以在 tspan 中指定对求解区间的分割,如:
[x,y]=ode23(fun,[0:0.1:0.5],1); %此时 x=[0:0.1:0.5]

matlab实例讲解欧拉法求解微分方程

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中使用欧拉法求解微分方程是非常简单而直观的。

Matlab求解微分方程

Matlab求解微分方程

y12
)
y2

y1
y1(0) 2, y2 (0) 0
1、建立m-文件vdp1000.m如下:
function dy=vdp1000(t,y)
dy=zeros(2,1);
dy(1)=y(2);
2 1.5
1 0.5
0 -0.5
-1 -1.5
-2 -2.5
0
500
1000
1500
2000
2500
3、使用泰勒公式
以此方法为基础,有龙格-库塔法、线性多步法等方 法。
4、数值公式的精度 当一个数值公式的截断误差可表示为O(hk+1)时
(k为正整数,h为步长),称它是一个k阶公式。 k越大,则数值公式的精度越高。
•欧拉法是一阶公式,改进的欧拉法是二阶公式。 •龙格-库塔法有二阶公式和四阶公式。 •线性多步法有四阶阿达姆斯外插公式和内插公式。
2、使用Matlab软件求数值解时,高阶微分方程必须 等价地变换成一阶微分方程组.
例4

d2x dt 2
1000(1
x
2
)
dx dt

x

0
x(0) 2; x'(0) 0
解: 令 y1=x,y2=y1’
则微分方程变为一阶微分方程组:

y1' y2

y2
'

1000(1
返回
微分方程的解析解
求微分方程(组)的解析解命令: dsolve(‘方程1’, ‘方程2’,…‘方程n’, ‘初始条件’, ‘自变量’)
记号: 在表达微分方程时,用字母 D 表示求微分,D2、D3 等

MATLAB数值分析实验五(欧拉法,荣格-库塔法解常微分方程)

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时,误差就很小了,再缩小一半,误差就更小了。

MATLAB解微分方程

MATLAB解微分方程
-2 -2.5
0
500
1000 1500 2000 2500 3000
数学建模实验项目2(9)
五、实验题目
实验要求:针对下列问题,先建立数学模型,再确定模型 中的参数,并有结论。
1、求下列微分方程的解析解
(1) y′ − 2 y = x2
(2) y′ − xy = x, y(0) = 1
(3) y′′ + y′ = x, y(0) = 1, y′(0) = 1
x2
)
dx dt

x
=
0
x(0) = 2; x'(0) = 0
解: 令 y1=x,y2=y1’
则微分方程变为一阶微分方程组:
y2
'
=
y1'= 1000(1−
y2 y12
)y2

y1
y1(0) = 2, y2(0) = 0
(1)建立m-文件vdp1000.m如下: function dy=vdp1000(t,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=1000*(1-y(1)^2)*y(2)-y(1);
其中k1是药物量被吸收到血液中的速率系数,k是血液中向体外 排除的速率系数,D是刚开始胃中或肌肉中的药物总量。 试用欧拉公式求上述微分方程数值解,并画出图形。(设
k1=0.6,k=0.2,D=200)
用Matlab软件求解微分方程 的解析解和数值解
数学与信息科学学院 孔祥庆
数学建模实验项目2(1)
一、实验名称:用Matlab软件求解微分方程的解析解和数值解 二、实验目的:
掌握用Matlab软件求解微分方程模型的解析解和数值解的方法 三、实验内容

用MATLAB函数编写并求解微分方程

用MATLAB函数编写并求解微分方程

机方针本质上解释微分方程的数值求解问题。
二、实验步骤
1. 建立数学模型,利用电路的拓扑元件的属性,列出网孔方程或节点方程,并简化为
标准形式的计算机可求解的一组微分方程组的过程。
2. 选择适合的计算机求解方法求解仿真模型。
3. 编写 MATLAB 仿真程序式建立 Simulink 模块方块图,调试并运行程序。
1.8 1.6 1.4 1.2
1 0.8 0.6 0.4 0.2
0 -0.2
-2
uC(t)
0
2
4
6
time sec
8
10
0.3 0.25
0.2 0.150.Fra bibliotek 0.050 -0.05
-0.1 -0.15
-0.2 -2
iL(t)
0
2
4
6
time sec
1.求解程序:零状态响应 %filename ex123.m L=1;c=0.1;for R=[1.5 3 5] [t,x]=ode 45(‘funcforex123’,[-1,10],[0;0],[],R,L,C); Figure(1);plot(t,x(:,1));hold onl xlibel(‘time sec’); Text(0.9,0.17,’\leftarrow i_L(t)’);grid; Figure(2);plot(t,x(:,2));hold on; xlabel(‘time sec’); Text(0.5,0.3,’\leftarrow u_C(t)’);grid;
用matlab函数编写并求解微分方程一实验原理为了对连续系统进行方针首先需要建立其数学模型然后利用计算机求这些数学模型从而得出数学模型的数值解

matlab欧拉折线法

matlab欧拉折线法

matlab欧拉折线法 欧拉折线法是一种常用的数值解法,用于求解微分方程的近似解。

在Matlab中,可以通过编写代码来实现欧拉折线法。

欧拉折线法的基本思想是将微分方程转化为差分方程,通过逐步逼近来计算函数的近似值。

具体步骤如下:
1. 确定微分方程的初值条件,即在某一点上的函数值和导数值。

2. 将自变量范围划分成若干个小区间,确定步长h。

3. 使用差分公式:y(i+1) = y(i) + h * f(x(i),
y(i)) 来计算下一个点的函数值,其中f(x,y)表示微分方程的右端项。

4. 重复步骤3,直到达到所需的自变量范围。

matlab 数学实验 实验报告 欧拉公式 ROSSLER微分方程

matlab 数学实验 实验报告 欧拉公式 ROSSLER微分方程

数学实验—实验报告一、实验项目: 二、实验目的和要求1、本章将对人口变化、动物种群变迁、网络系统的可靠性分析,介绍微分方程(组)的模型建立、数值解和图形解等方法,并用MATLAB 几何直观地展示各种求解方法的求解结果。

2、利用欧拉公式求解方程三、实验题目问题一:求微分方程的解析解,并画出它们的图形, y ’=y +2x , y (0)=1, 0<x <1;问题二:用向前欧拉公式和改进的欧拉公式求方程y ’=y -2x /y , y (0)=1的数值解(0≤x ≤1 , h =0.1) 要求编写程序。

问题三:Rossler 微分方程组当固定参数b=2, c=4时,试讨论随参数a 由小到大变化(如a ∈(0,0.65])而方程解的变化情况,并且画出空间曲线图形,观察空间曲线是否形成混沌状?问题四:水的流出时间一横截面积为常数A ,高为H 的水池内盛满水,由池底一横截面积为B 的小孔放水。

设水从小孔流出的速度为v=(2gh)0.5,求在任意时刻的水面高度和将水放空所需的时间。

时间t →高度h 。

问题五:考虑相互竞争模型两种相似的群体之间为了争夺有限的同一种事物来源和生存空间而进行生死存亡竞争时,往往是竞争力较弱的种群灭亡,而竞争力较强的种群达到环境容许的最大数量假设有甲、乙两个生物种群,当它们各自生存于一个自然环境中,均服从Logistics 规律。

三、实验过程问题一:用matlab 编写代码: x=[0,1]y=dsolve('Dy=y+2*x')y=dsolve('Dy=y+2*x', 'y(0)=1', 'x')ezplot(x,y)输出:y =-2*x+exp(t)*C1 (通解)y =-2*x-2+3*exp(x) 画图:x=0:0.01:1;y =-2*x-2+3*exp(x);plot(x,y)'''()x y z y x ayz b z x c =--⎧⎪=+⎨⎪=+-⎩0.10.20.30.40.50.60.70.80.91问题二: 1、分析:解:(1)解析解法得到其精确解:(2) 向前欧拉法:1(2/)n n n n n y y h y x y +=+- (1)2/n n nh y h x y =+- 迭代公式为 n+1y 1.10.2/n n n y x y =-,其中0y =y(0)=1(3)改进欧拉法:n+1nn n n n+1n +1n n n n n n n n n n n nn2n nnnn y =y +(h /2)*[(y -2x /y )+(y -2x /y )]=y +(h /2)*[(y -2x/y )+(y +h -2(x+h )/(yh ))] =y+(h /2)*[2y2x /y2(x +h )/(y+h )]=(1+h )y /2x/y(x +h )/(y+h )h h h h +--+-- 迭代公式为 n+1y 1.10.1/0.1()/()0.005n n n n n y x y x h y h =--+++,其中0y =y(0)=12、Matlab 编码x1(1)=0;y1(1)=1;y2(1)=1;h=0.1; for k=1:10 x1(k+1)=x1(k)+h;y1(k+1)=(1-h)*y1(k)+2*h*x1(k)/y1(k);y2(k+1)=(1+h)*y2(k)+(h*h)/2-h*x1(k)/y2(k)-h*(x1(k)+h)/(y2(k)+h); end x=0:0.1:1; y=(2*x+1).^(1/2);x1=x1(1:11),y=y(1:11),y1=y1(1:11),y2=y2(1:11), plot(x,y,x1,y1,'k:',x1,y2,'r--')显示图像及结果:x1 = 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000y = 1.0000 1.0954 1.1832 1.2649 1.3416 1.4142 1.4832 1.5492 1.6125 1.6733 1.7321y1 = 1.0000 0.9000 0.8322 0.7971 0.7926 0.8143 0.8557 0.9103 0.9731 1.0402 1.1092y2 = 1.0000 1.0959 1.1847 1.2679 1.3468 1.4222 1.4948 1.5653 1.6340 1.7016 1.768300.10.20.30.40.50.60.70.80.910.811.21.41.61.82图中,蓝色曲线是精确解,红色曲线是向前欧拉法曲线,黑色曲线是改进后欧拉法曲线问题三1、matlab 编程function r=rossler(t,x) global a; global b; global c;r=[-x(2)-x(3);x(1)+a*x(2);b+x(3)*(x(1)-c)];global a; global b; global c; b=2; c=4; t0=[0,200]; for a=0:0.02:0.65[t,x]=ode45('rossler',t0,[0,0,0]); subplot(1,2,1);plot(t,x(:,1),'r',t,x(:,2),'g',t,x(:,3),'b');title('x(ºìÉ«),y(ÂÌÉ«),z(ÀºÉ«)Ëæt±ä»¯Çé¿ö');xlabel('t'); subplot(1,2,2);plot3(x(:,1),x(:,2),x(:,3))title('Ïàͼ');xlabel('x');ylabel('y');zlabel('z'); pause end当a=0时,图像如下50100150200-1-0.8-0.6-0.4-0.200.20.40.6x(红色),y(绿色),z(篮色)随t 变化情况t-0.5相图z当a=0时,(x,y,z)收敛于(0,0.5,0.5)当a=0.12时,图像如下50100150200-1.2-1-0.8-0.6-0.4-0.200.20.40.6x(红色),y(绿色),z(篮色)随t 变化情况t-0.5相图z当a=0.28时,(x,y ,z)仍然收敛,但是收敛速度大大降低。

Matlab实验报告五(微分方程求解Euler折线法)

Matlab实验报告五(微分方程求解Euler折线法)
y=1;
szj=[x,y];
fori=1:n-1
y=y+h*subs(f,{'x','y'},{x,y});
x=x+h;
szj=[szj;x,y];
end
Szj;
plot(szj(:,1),szj(:,2))
3.结果
4.结论及分析
通过实验,结论正确,证明分析无误。
三、实验小结
这次实验说实话不是很难,因为这两个题型老师在课堂上已经讲得很清楚了,再说课本上还有类试题型,所以很轻松地做出来了。
y =
2^(1/2)/(4*exp(2^(1/2)*t)) - (2^(1/2)*exp(2^(1/2)*t))/4
x =
cosh(2^(1/2)*t) - (2^(1/2)*sinh(2^(1/2)*t))/2
y =
-(2^(1/2)*sinh(2^(1/2)*t))/2
4.结论及分析
通过实验,结论正确,证明分析无误。
问题2.用Euler折现法求解常微分方程 的数值解(步长 ),求解范围 ,并作出去图像.
1.分析问题
本题是用Euler折线法根据已知条件求解微分方程组的数值解,并作出它的图形。
2.问题求解
clear
f=sym('y-12*x^2/y^3');
a=0;
b=2;
h=0.001;
n=(b-a)/h+1;
x=0;
x=simple(x)
y=simple(y)
ezplot(x,y,[0,0.5]);axisauto
3.结果
x =
exp(2^(1/2)*t)/2 + 1/(2*exp(2^(1/2)*t)) - (2^(1/2)*exp(2^(1/2)*t))/4 + 2^(1/2)/(4*exp(2^(1/2)*t))

Matlab常微分方程的求解实验报告

Matlab常微分方程的求解实验报告

《数学实验》报告实验名称Mat lab常微分方程的求解学院计算机与通信工程学院专业班级计1103姓名------------- 学号-----------月6年2013.一、【实验目的】通过练习,熟悉Mat丄ab的求解常微分方程,函数文件的创建等。

了解Mat lab的命令窗口及其基本操作和常用命令。

通过练习,熟悉Mat lab的•些基本操作,掌握符号解法和数值解法,以及其中常用的方法。

二、【实验任务】1 > 求解微分力程y* =xsin (x) /cos (y)<»2、用数值方法求解下列微分方程,用不同颜色和线形将y和y,画在同•个图形窗口里:y f »+ty f-y=l-2t,初始时间:t=0;终止时间:t=n •初始条件:y |=0.1, y1|=0.2o €=oft=oo三、【实验程序】题一:y= dsolve('Dy=x*sin(x)/cos (y) ', 'x1)题二令:a=y, b二y,二a,, b'二y,,贝9 :a' =0*a+l*b+0;b,=l*a-t*b+(l-2*t);[a ;b,] = [0 1;1 —t][a;b] + [0;l](l-2*t);故化为一阶微分方程有:x =[0 1;1 -t]x+[0;l]u;其中:x'二[a' ;b'],u二(l~2*t) o初始条件:当t=0时,a=0. 1; b=0. 2oxd = mainfun( )function %UNTITLED Summary of this function goes here %ailed explanation goes here u=l-2 ・傑匸;xa=[0 1;1 -t][0;1]*u;endelf; tO=O; tf=pi;初始条fT%xOt=[0.1;0.2]; ,[t0r tf]/xOt) 'mainfun*[t,x]=ode23( y=x(:,1);Dy=x (:,2); ) ' 'r t, Dy, plot (t z y, ) 1 1Y 轴,匸轴,),ylabel (xlabel () 一阶导数'Dy''原始函数y \ legend (四、【实验结果】y =asin(C3 + sin(x) - x*cos(x))»funt =0.01450.08730.20140.32590.46210.61210.77780.96211.14821.27671.40531.51881.67061.86012.08912.35692.65462.96873.14160.1000 0.20000.1030 0.21580.1214 0.28830.159S 037980.2116 0.44790.2756 0.48470.3485 0.48130.4245 0.42820.4939 0.31680.5388 0.15990.5513 0.03200.5466 -0.10690.5272 -0.23580.4780 -0.41270.3788 -0.63470.2037 -0.8952-0.0742 ・ 1.1797-0.4677 -1.4646-0.9691 -1.7290-1.2793 -1.8586五、【实验总结】4通过本次实验,熟悉了M&t丄&b的一些基本操作,通过练习,会求解常微分方程,并熟练掌握符号解法和数值解法,以及其中常用的方法,包括欧拉方法、梯形公式、龙格库塔方法等。

微分方程的matlab求解

微分方程的matlab求解
2102111???????nyxfyxfhyynnnnnn改进欧拉公式??hhynhfxnyn主页下一页上一页??????????????21121211hkyxfkyxfkkkyynnnnnn返回数学实验x0精确解11004810187104081040810703110651148811966124931306613679向前欧拉111011029102910561109051131411783123051287413487向后欧拉改进欧拉11009110264105131051310830112091164512132126651324113855101020303040506070809110051019104121041210708110711149411972125001307213685步长h01的数值解比较表使用改进欧拉公式的例主页下一页上一页数学实验2龙格库塔法龙格库塔法是利用泰勒展式将yxh在x处展开并取其前面若干项来近似yxh而得到公式yxh?yxh??xyxh如果yxn?yn则yxn1的近似值为
数学实验
2、龙格-库塔法
龙格-库塔法是利用泰勒展式将y(x+h)在x处展开,并 取其前面若干项来近似y(x+h)而得到公式 y(x+h) y(x) + h j (x, y(x), h) 如果y(xn) yn,则y(xn+1)的近似值为:
yn+1 = yn + h j (xn, yn, h), n = 0, 1,…
精确解 1 1.0048
1.0187 1.0408 1.0703 1.1065 1.1488 1.1966 1.2493 1.3066 1.3679
上一页 下一页
向前欧拉 向后欧拉 1 1 1 1.0091

解微分方程欧拉法,R-K法及其MATLAB实例

解微分方程欧拉法,R-K法及其MATLAB实例

解微分方程的欧拉法,龙格-库塔法及其MATLAB简单实例欧拉方法(Euler method)用以对给定初值的常微分方程(即初值问题)求解分为前进EULER法、后退EULER法、改进的EULER法。

缺点:欧拉法简单地取切线的端点作为下一步的起点进行计算,当步数增多时,误差会因积累而越来越大。

因此欧拉格式一般不用于实际计算。

改进欧拉格式:为提高精度,需要在欧拉格式的基础上进行改进。

采用区间两端的斜率的平均值作为直线方程的斜率。

改进欧拉法的精度为二阶。

算法为:微分方程的本质特征是方程中含有导数项,数值解法的第一步就是设法消除其导数值。

对于常微分方程:x∈[a,b]y(a) = y0可以将区间[a,b]分成n段,那么方程在第xi点有y'(xi) = f(xi,y(xi)),再用向前差商近似代替导数则为:在这里,h是步长,即相邻两个结点间的距离。

因此可以根据xi点和yi点的数值计算出yi+1来:i=0,1,2,L这就是向前欧拉格式。

改进的欧拉公式:将向前欧拉公式中的导数f(xi,yi)改为微元两端导数的平均,即上式便是梯形的欧拉公式。

可见,上式是隐式格式,需要迭代求解。

为了便于求解,使用改进的欧拉公式:数值分析中,龙格-库塔法(Runge-Kutta)是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。

实际上,龙格-库塔法是欧拉方法的一种推广,向前欧拉公式将导数项简单取为f(xn,yn),而改进的欧拉公式将导数项取为两端导数的平均。

龙格-库塔方法的基本思想:在区间[xn,xn+1]内多取几个点,将他们的斜率加权平均,作为导数的近似。

龙格库塔法的家族中的一个成员如此常用,以至于经常被称为“RK4”或者就是“龙格库塔法”。

令初值问题表述如下。

则,对于该问题的RK4由如下方程给出:其中这样,下一个值(y n+1)由现在的值(y n)加上时间间隔(h)和一个估算的斜率的乘积决定。

该斜率是以下斜率的加权平均:k1是时间段开始时的斜率;k2是时间段中点的斜率,通过欧拉法采用斜率k1来决定y在点tn + h/2的值;k3也是中点的斜率,但是这次采用斜率k2决定y值;k4是时间段终点的斜率,其y值用k3决定。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

数学与信息科学系实验报告
实验名称微分方程求解
所属课程数学软件与实验
实验类型综合型实验
专业信息与计算科学
班级
学号
姓名
指导教师
一、实验概述【实验目的】 熟悉在Matlab 环境下求解常微分方程组和偏微分方程组的方法,掌握利用Matlab 软件进行常微分方程组和偏微分方程组的求解。

【实验原理】 1.dsolve(‘equ1’,’equ2’,...):matlab 求微分方程的解析解。

2.simplify(s):对表达式S 使用MAPLE 的化简规则进行化简。

3.[x,y]=dslove(‘方程1’,‘方程2’,...‘初始条件1’‘初始条件2’,..’自变量’):用字符串方程表示,自变量缺省值为t.4.ezplot(x,y,[tmin,tmax]):符号函数的作图命令。

【实验环境】 MatlabR2010b 二、实验内容问题1. 求微分方程组在初始条件下的解,并00dx x y dt dy x y dt ⎧++=⎪⎪⎨⎪+-=⎪⎩00|1,|0t t x y ====[0,0.5]t ∈画出函数的图像. ()y f x =1.分析问题本题是根据初始条件求微分方程组的特解,并根据t 的范围画出函数的图形。

2.问题求解
syms x y t [x,y]=dsolve('Dx+x+y=0','Dy+x-y=0','x(0)=1','y(0)=0','t')x=simple(x)y=simple(y)ezplot(x,y,[0,0.5]);axis auto 3.结果x =exp(2^(1/2)*t)/2 + 1/(2*exp(2^(1/2)*t)) -
(2^(1/2)*exp(2^(1/2)*t))/4 + 2^(1/2)/(4*exp(2^(1/2)*t))
y =2^(1/2)/(4*exp(2^(1/2)*t)) - (2^(1/2)*exp(2^(1/2)*t))/4x =cosh(2^(1/2)*t) - (2^(1/2)*sinh(2^(1/2)*t))/2。

相关文档
最新文档