Matlab 解常微分方程的初值问题
重要：MATLAB常微分方程(组)数值解法
Matlab常微分方程求解问题分类
边值问题：
初值问题：
• 定解附加条件在自变量 的一端
• 一般形式为： y' f (x, y)
y(a)
y0
• 初值问题的数值解法一 般采用步进法，如 Runge-Kutta法
➢ 在自变量两端均给定附加 条件
y' f (x, y)
➢ 一般形式：y(a)y1, y(b)y2
1.根据常微分方程要求的求解精度与速度要求
求解初值问题：
y
'
y
2x y
y ( 0 ) 1
(0x1)
比较ode45和ode23的求解精度和速度
ode45和ode23的比较-1
function xODE clear all clc
format long
y0 = 1; [x1,y1] = ode45(@f,[0,1],y0); [x2,y2] = ode23(@f,[0,1],y0); plot(x1,y1,'k-',x2,y2,'b--') xlabel('x') ylabel('y')
rD = k(3)*C(2)-k(5)*C(4);
rE = k(4)*C(3)+k(5)*C(4);
% Mass balances dCdt = [rA; rB; rC; rD; rE];
三个串联的CSTR等温反应器（例4-3）
function IsothermCSTRs clear all clc CA0 = 1.8; % kmol/m^3 CA10 = 0.4; % kmol/m^3 CA20 = 0.2; % kmol/m^3 CA30 = 0.1; % kmol/m^3 k = 0.5; % 1/min tau = 2; stoptime = 2.9; % min [t,y] = ode45(@Equations,[0 stoptime],[CA10 CA20 CA30],[],k,CA0,tau); disp(' Results:') disp(' t CA1 CA2 CA3') disp([t,y]) plot(t,y(:,1),'k--',t,y(:,2),'b:',t,y(:,3),'r-') legend('CA_1','CA_2','CA_3') xlabel('Time (min)') ylabel('Concentration') % -----------------------------------------------------------------function dydt = Equations(t,y,k,CA0,tau) CA1 = y(1); CA2 = y(2); CA3 = y(3); dCA1dt = (CA0-CA1)/tau - k*CA1; dCA2dt = (CA1-CA2)/tau - k*CA2; dCA3dt = (CA2-CA3)/tau - k*CA3; dydt = [dCA1dt; dCA2dt; dCA3dt];
四阶龙格-库塔法求解常微分方程的初值问题-matlab通用程序
参考教材《数值分析》李乃成.梅立泉clearclcformat longm=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);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)); %计算K1x=x+h/2;y(1,1)=y1+K(1,1)/2;y1=y(1,1);K(1,2)=h*(eval(ym)); %计算K2x=x;y(1,1)=y1+K(1,2)/2-K(1,1)/2;y1=y(1,1);K(1,3)=h*(eval(ym)); %计算K3x=x+h/2;y(1,1)=y1+K(1,3)-K(1,2)/2;y1=y(1,1);K(1,4)=h*(eval(ym)); %计算K4y11(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;endy11else %高阶初值问题mm=(b-a)/h; %一共要求解mm个数据点for k2=1:m %读取初值条件fprintf('请输入%d阶导数的初值f(%d)(a)=\n',(k2-1),(k2-1));y(k2,1)=input('=');endfor k2=1:my22(1,k2)=y(k2,1); %先把初值保存在矩阵y22(m,n)中,m表示第几个所求点,n表示第n阶初值endx=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表示K1endK(m,1)=h*(eval(ym));x=x+h/2; %求解K1之前,先重新对x和y赋值for k3=1:my(k3,1)=y(k3,1)+K(k3,1)/2;endfor k=1:(m-1) %计算K2K(k,2)=h*y(k+1,1);endK(m,2)=h*(eval(ym));x=x;for k3=1:my(k3,1)=y(k3,1)-K(k3,1)/2+K(k3,2)/2;endfor k=1:(m-1) %计算K3K(k,3)=h*y(k+1,1);endK(m,3)=h*(eval(ym));x=x+h/2;for k3=1:my(k3,1)=y(k3,1)+K(k3,3)-K(k3,2)/2; %这里容易出错endfor k=1:(m-1) %计算K4K(k,4)=h*y(k+1,1);endK(m,4)=h*(eval(ym));for k5=1:my22(k4,k5)=y22(k4-1,k5)+(K(k5,1)+2*K(k5,2)+2*K(k5,3)+K(k5,4))/6; %这里,除了要求出下一个点的数值,还要求出相应的导数值endfor k6=1:m %除了对y(1,1)重新赋值外,还要对y(2,1)等重新赋值y(k6,1)=y22(k4,k6);endx=a+(k4-1)*h;endy22(:,1) end。
radau伪谱法 matlab
Radau伪谱法是一种在Matlab中常用的数值方法，用于解决特定类型的常微分方程初值问题。
这种方法通过将微分方程离散化为代数方程组，并使用Radau伪谱项目进行数值求解。
下面将介绍Radau伪谱法的基本原理、Matlab中的实现以及其在工程实践中的应用。
一、Radau伪谱法的基本原理1.1 Radau伪谱法的概念Radau伪谱法是一种基于伪谱分解的高阶数值方法，用于求解常微分方程的初值问题。
它在时间上采用隐式Runge-Kutta方法，并利用伪谱投影将微分方程离散化为代数方程组，进而通过求解代数方程组来获得微分方程的数值解。
1.2 Radau伪谱法的数学原理Radau伪谱法的数学原理涉及到伪谱变换和伪谱求解两个基本步骤。
首先通过伪谱变换将微分方程转化为代数方程组，然后通过求解代数方程组来获得微分方程的数值解。
1.3 Radau伪谱法的优点与传统的数值方法相比，Radau伪谱法具有更高的收敛速度和更好的数值稳定性。
特别是对于一些刚性和高阶微分方程问题，Radau伪谱法表现出了更好的数值效果。
二、Matlab中Radau伪谱法的实现2.1 Matlab中的Radau伪谱工具箱Matlab提供了丰富的数值计算工具箱，其中包括了Radau伪谱工具箱，用户可以直接调用该工具箱来实现Radau伪谱法的数值计算。
2.2 Radau伪谱法在Matlab中的调用方法在Matlab中，用户可以通过调用Radau伪谱工具箱中的相关函数来进行Radau伪谱法的数值计算。
用户需要提供待求解的微分方程、初值条件以及数值求解的参数设置，然后调用相应的函数即可获得微分方程的数值解。
2.3 Matlab中Radau伪谱法的优化和加速Matlab中的Radau伪谱工具箱经过了多次优化和加速，能够更加高效地求解微分方程初值问题。
用户可以通过调整参数设置和算法选项来获得更高质量的数值解，并在工程实践中取得更好的效果。
三、Radau伪谱法在工程实践中的应用3.1 Radau伪谱法在航天工程中的应用在航天工程中，常常需要对复杂的动力学系统进行数值仿真和数值控制设计。
matlab_常微分方程数值解法
dt 2
简朴问题可以求得解析解，多数实际问题靠数值求解 。
第4页
一阶常微分方程(ODE )初值问题 ： ODE ：Ordinary Differential Equation
dy
f
(x,
y)
dx
x0 x xn
y(x0 ) y0
数值解法就是求y(x)在某些分立旳节点 xn 上旳近似值 yn，用以近似y(xn)
x0
y0
x1 f y(x), x dx
x0
x2 f y(x), x dx
x1
y(x1) f y(x1), x1 h
第17页
同样，在[x0,xn+1] ，积分采用矩形近似，得：
y(xn1) y0
f xn1
x0
y(x), x dx
y(xn ) f y(xn ), xn h
yn y(xn )
第5页
2、欧拉近似办法
2.1 简朴欧拉(L.Euler, 1707-1783)办法。
dy
dx
f
(y, x)
y(x0 ) y0
欧拉数值算法就是由初值通过递推求解，递推求解
就是从初值开始，后一种函数值由前一种函数值得到。核 心是构造递推公式。
y0 y1 y2 yn
第6页
i 1,2,...
第36页
没有一种算法可以有效地解决所有旳 ODE 问题，因此 MATLAB 提供了多种ODE函数。
函数 ODE类
特点
阐明
型
ode45
非刚性 单步法；4，5 阶 R-K 措施；合计 大部分场合旳首选措施
截断误差为 (△x)3
ode23
非刚性 单步法；2，3 阶 R-K 措施；合计 使用于精度较低旳情形
MATLAB常微分方程数值解——欧拉法、改进的欧拉法与四阶龙格库塔方法
MATLAB常微分⽅程数值解——欧拉法、改进的欧拉法与四阶龙格库塔⽅法MATLAB常微分⽅程数值解作者：凯鲁嘎吉 - 博客园1.⼀阶常微分⽅程初值问题2.欧拉法3.改进的欧拉法4.四阶龙格库塔⽅法5.例题⽤欧拉法，改进的欧拉法及4阶经典Runge-Kutta⽅法在不同步长下计算初值问题。
步长分别为0.2,0.4,1.0.matlab程序：function z=f(x,y)z=-y*(1+x*y);function R_K(h)%欧拉法y=1;fprintf('欧拉法：x=%f, y=%f\n',0,1);for i=1:1/hx=(i-1)*h;K=f(x,y);y=y+h*K;fprintf('欧拉法：x=%f, y=%f\n',x+h,y);endfprintf('\n');%改进的欧拉法y=1;fprintf('改进的欧拉法：x=%f, y=%f\n',0,1);for i=1:1/hx=(i-1)*h;K1=f(x,y);K2=f(x+h,y+h*K1);y=y+(h/2)*(K1+K2);fprintf('改进的欧拉法：x=%f, y=%f\n',x+h,y);endfprintf('\n');%龙格库塔⽅法y=1;fprintf('龙格库塔法：x=%f, y=%f\n',0,1);for i=1:1/hx=(i-1)*h;K1=f(x,y);K2=f(x+h/2,y+(h/2)*K1);K3=f(x+h/2,y+(h/2)*K2);K4=f(x+h,y+h*K3);y=y+(h/6)*(K1+2*K2+2*K3+K4);fprintf('龙格库塔法：x=%f, y=%f\n',x+h,y);end结果：>> R_K(0.2)欧拉法：x=0.000000, y=1.000000欧拉法：x=0.200000, y=0.800000欧拉法：x=0.400000, y=0.614400欧拉法：x=0.600000, y=0.461321欧拉法：x=0.800000, y=0.343519欧拉法：x=1.000000, y=0.255934改进的欧拉法：x=0.000000, y=1.000000改进的欧拉法：x=0.200000, y=0.807200改进的欧拉法：x=0.400000, y=0.636118改进的欧拉法：x=0.600000, y=0.495044改进的欧拉法：x=0.800000, y=0.383419改进的欧拉法：x=1.000000, y=0.296974龙格库塔法：x=0.000000, y=1.000000龙格库塔法：x=0.200000, y=0.804636龙格库塔法：x=0.400000, y=0.631465龙格库塔法：x=0.600000, y=0.489198龙格库塔法：x=0.800000, y=0.377225龙格库塔法：x=1.000000, y=0.291009>> R_K(0.4)欧拉法：x=0.000000, y=1.000000欧拉法：x=0.400000, y=0.600000欧拉法：x=0.800000, y=0.302400改进的欧拉法：x=0.000000, y=1.000000改进的欧拉法：x=0.400000, y=0.651200改进的欧拉法：x=0.800000, y=0.405782龙格库塔法：x=0.000000, y=1.000000龙格库塔法：x=0.400000, y=0.631625龙格库塔法：x=0.800000, y=0.377556>> R_K(1)欧拉法：x=0.000000, y=1.000000欧拉法：x=1.000000, y=0.000000改进的欧拉法：x=0.000000, y=1.000000改进的欧拉法：x=1.000000, y=0.500000龙格库塔法：x=0.000000, y=1.000000龙格库塔法：x=1.000000, y=0.303395注意：在步长h为0.4时，要将for i=1:1/h改为for i=1:0.8/h。
实验报告七常微分方程初值问题的数值解法
浙江大学城市学院实验报告课程名称数值计算方法实验项目名称常微分方程初值问题的数值解法 实验成绩指导老师签名日期2015/12/16 一.实验目的和要求1． 用Matlab 软件掌握求微分方程数值解的欧拉方法和龙格－库塔方法； 2． 通过实例学习用微分方程模型解决简化的实际问题;二.实验内容和原理编程题2-1要求写出Matlab 源程序m 文件,并有适当的注释语句；分析应用题2-2,2-3,2-4,2-5要求将问题的分析过程、Matlab 源程序和运行结果和结果的解释、算法的分析写在实验报告上; 2-1 编程编写用向前欧拉公式和改进欧拉公式求微分方程数值解的Matlab 程序,问题如下：在区间[],a b 内(1)N +个等距点处,逼近下列初值问题的解,并对程序的每一句添上注释语句; Euler 法y=eulera,b,n,y0,f,f1,b1改进Euler 法y=eulerproa,b,n,y0,f,f1,b1 2-2 分析应用题假设等分区间数100n =,用欧拉法和改进欧拉法在区间[0,10]t ∈内求解初值问题()()20(0)10y t y t y '=-⎧⎨=⎩并作出解的曲线图形,同时将方程的解析解也画在同一张图上,并作比较,分析这两种方法的精度; 2-3 分析应用题用以下三种不同的方法求下述微分方程的数值解,取10h = 画出解的图形,与精确值比较并进行分析; 1欧拉法； 2改进欧拉法； 3龙格－库塔方法；2-4 分析应用题考虑一个涉及到社会上与众不同的人的繁衍问题模型;假设在时刻t 单位为年,社会上有人口()x t 人,又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人;而固定比例为r 的所有其他的后代也是与众不同的人;如果对所有人来说出生率假定为常数b ,又如果普通的人和与众不同的人的婚配是任意的,则此问题可以用微分方程表示为：其中变量()()()i p t x t x t =表示在时刻t 社会上与众不同的人的比例,()i x t 表示在时刻t 人口中与众不同的人的数量;1假定(0)0.01,0.02p b ==和0.1r =,当步长为1h =年时,求从0t =到50t =解()p t 的近似值,并作出近似解的曲线图形;2精确求出微分方程的解()p t ,并将你当50t =时在分题b 中得到的结果与此时的精确值进行比较; MATLAB 相关函数求微分方程的解析解及其数值的代入dsolve‘egn1’,‘egn2’,‘x ’ subsexpr,{x,y,…},{x1,y1,…}其中‘egn i ’表示第i 个方程,‘x ’表示微分方程中的自变量,默认时自变量为t ; subs 命令中的expr 、x 、y 为符合型表达式,x 、y 分别用数值x1、x2代入; >>symsxyz>>subs'x+y+z',{x,y,z},{1,2,3} ans= 6>>symsx>>subs'x^2',x,2 ans= 4>>s=dsolve‘12Dy y ∧=+’,‘(0)1y =’,‘x ’ ans= >>symsx >>subss,x,2 ans=右端函数(,)f x y 的自动生成f=inline ‘expr ’,’var1’,‘var2’,……其中’expr ’表示函数的表达式,’var1’,‘var2’表示函数表达式中的变量,运行该函数,生成一个新的函数表达式为fvar1,var2,……; >>f=inline'x+3y','x','y' f=Inlinefunction: fx,y=x+3y >>f2,3 ans= 114,5阶龙格－库塔方法求解微分方程数值解t,x=ode45f,ts,x0,options其中f 是由待解方程写成的m 文件名；x0为函数的初值；t,x 分别为输出的自变量和函数值列向量,t的步长是程序根据误差限自动选定的;若ts=t0,t1,t2,…,tf,则输出在自变量指定值,等步长时用ts=t0:k:tf,输出在等分点；options 用于设定误差限可以缺省,缺省时设定为相对误差310-,绝对误差610-,程序为：options=odeset ‘reltol ’,rt,’abstol ’,at,这里rt,at 分别为设定的相对误差和绝对误差;常用选项见下表;选项名 功能 可选值 省缺值 AbsTol 设定绝对误差正数 RelTol 设定相对误差 正数InitialStep 设定初始步长 正数 自动 MaxStep设定步长上界正数MaxOrder 设定ode15s 的最高阶数 1,2,3,4,5 5 Stats 显示计算成本统计 on,off off BDF 设定ode15s 是否用反向差分on,offoff例：在命令窗口执行>>odefun =inline ‘2*y t y -’,‘t ’,‘y ’；>>[],45(,[0,4],1)t y ode odefun =；ans=>>t y ‘o-’,%解函数图形表示>>45(,[0,4],1)ode odefun %不用输出变量,则直接输出图形 >>[],45(,0:4,1)t y ode odefun =；[],t yans=三.操作方法与实验步骤包括实验数据记录和处理2-1编程编写用向前欧拉公式和改进欧拉公式求微分方程数值解的Matlab 程序,问题如下：在区间[],a b 内(1)N +个等距点处,逼近下列初值问题的解,并对程序的每一句添上注释语句; Euler 法y=eulera,b,n,y0,f,f1,b1改进Euler 法y=eulerproa,b,n,y0,f,f1,b1Euler 法y=eulera,b,n,y0,f,f1,b1 y=zeros1,n+1; y1=y0; h=b-a/n; x=a:h:b; fori=1:n; yi+1=yi+hfxi,yi; end plotx,y holdon%求微分方程的精确解 x1=linspacea,b,100; '精确解为' s=dsolvef1,b1,'x' symsxy1=zeros1,100; for i=1:100y1i=subss,x,x1i; endplotx1,y1,'r'title'红色代表精确解'改进Euler 法y=eulerproa,b,n,y0,f,f1,b1 %求微分方程的数值解 y=zeros1,n+1; y1=y0; h=b-a/n; x=a:h:b; fori=1:n; T1=fxi,yi; T2=fxi+1,yi+hT1; yi+1=yi+h/2T1+T2; end plotx,y holdon%求微分方程的精确解 x1=linspacea,b,100; '精确解为' s=dsolvef1,b1,'x' symsxy1=zeros1,100; fori=1:100 y1i=subss,x,x1i; endplotx1,y1,'r'title'红色代表精确解' 2-2分析应用题假设等分区间数100n =,用欧拉法和改进欧拉法在区间[0,10]t ∈内求解初值问题()()20(0)10y t y t y '=-⎧⎨=⎩并作出解的曲线图形,同时将方程的解析解也画在同一张图上,并作比较,分析这两种方法的精度;1向前欧拉法>>euler0,10,100,10,inline'y-20','x','y','Dy=y-20','y0=10' ans= 精确解为 s= 20-10expx ans= +005Columns1through8（2）改进欧拉法>>eulerpro0,10,100,10,inline'y-20','x','y','Dy=y-20','y0=10' ans= 精确解为 s= 20-10expx ans= +005Columns1through8改进欧拉法的精度比向前欧拉法更高; 2-3分析应用题用以下三种不同的方法求下述微分方程的数值解,取10h = 画出解的图形,与精确值比较并进行分析; 1欧拉法； 2改进欧拉法；2-4分析应用题考虑一个涉及到社会上与众不同的人的繁衍问题模型;假设在时刻t 单位为年,社会上有人口()x t 人,又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人;而固定比例为r 的所有其他的后代也是与众不同的人;如果对所有人来说出生率假定为常数b ,又如果普通的人和与众不同的人的婚配是任意的,则此问题可以用微分方程表示为：其中变量()()()i p t x t x t =表示在时刻t 社会上与众不同的人的比例,()i x t 表示在时刻t 人口中与众不同的人的数量;1假定(0)0.01,0.02p b ==和0.1r =,当步长为1h =年时,求从0t =到50t =解()p t 的近似值,并作出近似解的曲线图形;2精确求出微分方程的解()p t ,并将你当50t =时在分题b 中得到的结果与此时的精确值进行比较;1>>euler0,50,50,,inline'','t','p','Dp=','p0= 1' ans= 精确解为 s=1-99/100expx/500 ans=Columns1through82>>dsolve'Dp=','p0=','t' ans=1-99/100expt/500 >>1-99/100exp ans=与欧拉法求得的精确值差0,0001四.实验结果与分析。
matlab解常微分方程
matlab解常微分⽅程1. ODE常微分⽅程ordinary differential equation的缩写，此种表述⽅式常见于编程，如MATLAB中Simulink求解器solver已能提供了7种微分⽅程求解⽅法：ode45(Dormand-Prince)，ode23(Bogacki-Shampine)，ode113(Adams)，ode15s(stiff/NDF)，ode23s(stiff/Mod. Rosenbrock)，ode23t(mod.stiff/Trapezoidal)，ode23tb(stiff/TR-BDF2)。
微分⽅程、微分⽅程组⾃标量 因变量 ⼀元 多元 函数 映射⼀元：只有⼀个因变量多元：有多个因变量导数 偏导：谁对谁的导数，因变量对⾃变量的导数，默认或缺省⾃变量为t 、x ？⼀元⽅程 多元⽅程 多元⽅程组 n个⽅程解n个未知量微分⽅程 ⼀阶 ⾼阶微分⽅程 ⼀阶微分⽅程组⼀阶常微分⽅程：Dx/dt + x = e^t⾼阶常微分⽅程：d^2x/dt^2+dx/dt+x=e^2t⼀阶微分⽅程组（多元）：dy/dt+x=e^2tdx/dt+2y-x=e^t初始条件：dy/dt0=... dx/dt0=... y0=... x0=...可以解出：y=f(t)=.... x=f(t)=.... 两个⽅程解两个未知数（因变量）⼀个N阶（多元）微分⽅程可以写成（分解成）N个⼀阶微分⽅程（即微分⽅程组）如：x.. + 2x. -x = u令x.=x2; x=x1 则...微分⽅程的精确解： r=dsolve('eqn1','eqn2',...,'cond1','cond2',...,'var').数值解： [t,y]=solver('odefun',tspan,y0,options)1. 求精确解1.微分⽅程r=dsolve('eqn1','eqn2',...,'cond1','cond2',...,'var').该命令中可以⽤D表⽰微分符号，其中D2表⽰⼆阶微分，D3表⽰三阶微分，以此类推。
带约束的常微分方程组matlab求解初值问题
带约束的常微分方程组matlab求解初
值问题
常微分方程组的初值问题是一种常见的数学问题，通过使
用Matlab可以方便地对其进行求解。
在求解常微分方程组时，如果存在约束条件，可以通过添加约束条件来求解。
需要定义常微分方程组的函数形式，以及约束条件的函数
形式。
然后，可以利用Matlab中的ode45函数来求解带约束
的常微分方程组的初值问题。
在使用ode45函数时，需要将约束条件作为一个输入参数
传入，以确保在求解过程中满足约束条件。
通过调整初始条件和约束条件，可以得到满足要求的数值解。
带约束的常微分方程组的初值问题是一种常见且重要的数
学问题，通过Matlab可以较为方便地进行求解。
在实际操作中，需要注意对约束条件的处理，以确保求解过程中满足约束条件，得到准确的数值解。
希望以上方法对您有所帮助。
如果还有其他问题，欢迎继续咨询。
感谢您的提问。
第5次实验报告(常微分方程初值问题的数值解法)
班级： 学号： 姓名： 成绩：实验5 常微分方程初值问题的数值解法实验1实验目的1）熟悉欧拉法、改进欧拉法和龙格-库塔法的原理。
2）根据以上方法，编程求解常微分方程初值问题。
2 实验内容（1）编写程序，用以上各种方法求解教材P232例7-1、习题6、11的初值问题。
(2) 使用系统自带的函数dsolve 和ode45求例7-1的符号解析解和数值解。
3实验原理求解微分方程初值问题00(,)()y f x y y x y '=⎧⎨=⎩ （1） 欧拉法（显式）：10(,)n n n n n y y hf x y x x nh +=+⎧⎨=+⎩（2） 改进欧拉法(0)1(0)111(,)(,)(,)2n n n n n n n n n n y y hf x y h y y f x y f x y ++++⎧=+⎪⎨⎡⎤=++⎪⎣⎦⎩ （3） 经典龙格-库塔法（四阶）11234121324300(22)6(,)(,)22(,)22(,)()i i i i i i i i i i h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK y y x +⎧=++++⎪⎪=⎪⎪=++⎪⎨⎪=++⎪⎪⎪=++⎪=⎪⎩4实验步骤1)建立函数文件，根据各公式编写程序；2)上机调试程序，运行程序进行计算，记录计算结果；3)分析各公式计算结果，比较各公式的优缺点。
5 程序设计欧拉法改进欧拉法function Euler1(x0,y0,h,n)%(x0,y0)：方程的初值%h：步长%n：计算的步数for i=1:nx=x0+h;yp=y0+h*f(x0,y0);yc=y0+h*f(x,yp);y=(yp+yc)/2;x %在屏幕显示每一步的x值y %在屏幕显示每一步计算的方程的数值解 x0=x;y0=y;end经典龙格-库塔法1）函数function f=f(x,y) f=y-2*x/y;end6总结注：若要更改matlab计算的数值类型，可以通过在matlab中设置实现：File -> Preferences ->Array Editor窗口中，Format 下方将Default array format设置为：long解微分方程的MATLAB命令MATLAB中主要用dsolve求微分方程的符号解析解，ode45求数值解。
Matlab微分方程的解法
-0.5
-0.55
-0.6
-0.65
-0.7
-0.75
-0.8
-0.85
-0.9
-0.95
-1
0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
1
time t0=0,tt=1
图3 给定新的初始数据，由函数xprim2定义的ODE解的图形
（d） 求解下面方程组并不难：
x x x x ì ' = - 0.1
在下面的初值问题中，有两个未知函数：x1(t)和x2(t),并用以下式子表达其微... 页码，1/11
Matlab关于微分方程的解法
MATLAB使用龙格-库塔-芬尔格（Runge-Kutta-Fehlberg）方法来解ODE问题。在有限点内计算求解。而 这些点的间距有解的本身来决定。当解比较平滑时，区间内使用的点数少一些，在解变化很快时，区间内应使 用较多的点。 为了得到更多的有关何时使用哪种解法和算法的信息，推荐使用helpdesk。所有求解方程通用的语法或句法在 命令集中头两行给出。时间间隔将以向量t=[t0,tt]给出。 命令ode23可以求解(2,3)阶的常微分方程组，函数ode45使用(4,5)阶的龙格-库塔-芬尔格方法。注意，在这种情 况下x’是x的微分不是x的转置。 在命令集中solver将被诸如ode45函数所取代。
0.6
0.55
0.5
0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
1
time t0=0,tt=1
图1 由函数xprim1定义的ODE解的图形
（b） 解下面的ODE过程是等价的：
ïíìx' = x2
ïîx(0) = 1
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
题目：Matlab 解常微分方程的初值问题
一、 设计目的：
1、熟练掌握Matlab 的基本编程方法，及其编程风格。
2、熟练掌握Matlab 常用函数的使用。
3、与本专业相关知识相结合，掌握其在程序开发中的应用方法 以及和word 、C 语言等接口方法。
4、通过计算机数值求解的方式来加深微分方程解的理解。
5、熟悉初等方法可获得解析解之外的数值近似解的求解方法，提 高对差分格式的认识和离散化分析问题的技巧，加深对理论课程的学习和理解，为数学专业和信息与计算科学专业其他后继课程的学习打好基础。
二、 设计内容：
已知一个三阶微分方程:012'
'
'2'
''=+⎪⎭⎫ ⎝⎛--y y y y
y ，利用matlab
软件求这个三阶微分方程在初值 ()()()1
0001
0'
''
-===y y y 下的解。
原三阶微分方程可化为：y y y y y '
'
'2
'
''12-⎪⎭
⎫ ⎝⎛-=
令
y y y y y
y '
''
321==
= 则原三阶微分方程可化为微分方程组213123
213212''
'
y y y y y y y y y -⎪⎭
⎫ ⎝⎛-=== 在初值 ()()()10001
0'
''
-===y y y 下的解。
三、 程序流程：
四、程序代码：
%编写函数文件rigid.m
function dy = rigid(t,y)
dy = zeros(3,1); % a column vector dy(1) = y(2) ;
dy(2) =y(3);
dy(3) = 2*(1-y(1)^2)*y(3)-y(1)*y(2);
%调用函数ode45求解，时间区间为[0,10]
[t,Y] = ode45(@rigid,[0 10],[1 0 -1])
t =
0.0001
0.0001
0.0002
0.0002
0.0005
0.0007
………………
0.9383
1.0665
1.1947
1.2918
1.3889
1.4860
………………
6.2916
6.2922
6.2928
6.2934
6.2940
6.2947
6.
由于数据太多，这里只列举部分
%绘制解的曲线
plot(t,Y(:,1),'-',t,Y(:,2),'-.',t,Y(:,3),'.')
%给图形加标注
title('Solution of Rigid Equation') xlabel('time T')
ylabel('solution Y')
legend('Y1','Y2','Y3')
012345678910
未加图形标注时的图
01234
5678910
Solution of Rigid Equation
time T
s o l u t i o n Y
加了图形标注后的图
输出结果[T,Y]中T为时间点组成的向量。
Y为对应于T中时间点的y(1)、y(2)和y(3)的值。
此次利用matlab数值方法来求解微分方程主要是把求解的时间划分成有限步，对应于每一步将计算出一个解，如果求得的解不满足误差限制，则减少步长，再求解。
如此重复，直到满足误差限为止。
七、课设总结：
最初拿到题目后，开始读题，知道题目所表达的意思及我们所要完成的目的，达到的效果后，便开始了做题。
首先是找一个三阶微分方程。
然后把它化简为标准形式，再利用matlab软件求解。
其中，在利用matlab求解时遇到一些问题，比如画图时调用已经编号的rigid函数时的调用格式不正确，还有就是给图形家标注时程序的引号没有切换成英文输入法状态下的等一些问题。
但是经过多次调试，检查，修改后程序运行总算成功。
通过此次课程设计，我已经基本掌握Matlab的基本编程方法，及其编程风格。
且能较熟练掌握Matlab常用函数的使用。
此次设计后我的收获不少。
八、参考文献：
[1] 张圣勤编 MATLAB7.0 机械工业出版社
[2]周义仓靳祯秦军林编常微分方程极其应用科学出版社
[3]韩明王家宝李林编数学实验（matlab版）同济大学出版社
[4]汪晓银皱庭荣编数学软件与数学实验科学出版社。