常微分方程的初值问题_Matlab代码
- 格式:pdf
- 大小:130.09 KB
- 文档页数:11
参考教材《数值分析》李乃成.梅立泉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。
常微分方程MATLAB程序以下是一个简单的MATLAB 程序,用于求解一阶常微分方程:matlab复制代码% 定义微分方程 dy/dx = f(x, y)f = @(x, y) -x*y;% 初始条件 y(0) = 1y0 = 1;% 定义 x 的范围xspan = [0, 10];% 使用 MATLAB 内置函数 ode45 进行求解[t, y] = ode45(f, xspan, y0);% 绘制解的图形plot(t, y(:,1));xlabel('x');ylabel('y');title('Solution of the differential equation dy/dx = -xy');在这个程序中,我们定义了一个一阶常微分方程dy/dx = -xy,并使用MATLAB 内置函数ode45进行求解。
初始条件为y(0) = 1,求解范围为xspan = [0, 10]。
最后,我们使用plot函数绘制了解的图形。
这个程序是用来求解一阶常微分方程的,而这个方程是dy/dx = -xy。
这是一个简单的线性方程,但它的解在物理和工程中有许多实际应用。
接下来,我们逐行解释一下代码:1.% 定义微分方程 dy/dx = f(x, y):这是一个注释,说明下面的代码是定义微分方程。
2. f = @(x, y) -x*y;:这行定义了一个匿名函数f,它接受两个参数x和y,并返回-x*y。
这个函数就是我们的微分方程dy/dx的右边部分。
3.% 初始条件 y(0) = 1:这是一个注释,说明下面的代码是定义初始条件。
4.y0 = 1;:这行定义了初始条件y(0) = 1,也就是说当x=0时,y=1。
5.% 定义 x 的范围:这是一个注释,说明下面的代码是定义自变量x的范围。
6.xspan = [0, 10];:这行定义了自变量x的范围从0到10。
7.% 使用 MATLAB 内置函数 ode45 进行求解:这是一个注释,说明下面的代码将使用MATLAB 的内置函数ode45来求解微分方程。
浙江大学城市学院实验报告课程名称数值计算方法实验项目名称常微分方程初值问题的数值解法 实验成绩指导老师签名日期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. 亚当斯(Adams)预测-修正算法由亚当斯-巴什福特(Adams-Bashforth)显式预测公式和亚当斯-莫顿(Adams-Moulton)隐式修正公式组成的预测-修正(PECE)对。
function [ YMat ] = Adams( func, tvec, y_init, order )% Adams预测-修正算法,⽤于求解常微分初值问题% 输⼊四个参数:函数句柄func(接收列向量、返回列向量),积分时间列向量tvec,初值⾏向量y_init,阶数order;% 输出⼀个参数:数值解,每⼀⾏对应积分时间列向量的⼀⾏,各列为变量⼀个分量。
switch ordercase '4'row = size(tvec, 1); col = size(y_init, 2);YMat = zeros(row, col);YMat(1:4, :) = Runge_Kutta(func, tvec(1:4), y_init, '4');for i=4:row - 1stepsize = tvec(i + 1) - tvec(i);ydiff0 = func(tvec(i), YMat(i, :).');ydiff1 = func(tvec(i - 1), YMat(i - 1, :).');ydiff2 = func(tvec(i - 2), YMat(i - 2, :).');ydiff3 = func(tvec(i - 3), YMat(i - 3, :).');y_predict = YMat(i, :).' + (55*ydiff0 - 59*ydiff1 + 37*ydiff2 - 9*ydiff3)*stepsize/24;y_corrector = YMat(i, :).' + (9*func(tvec(i + 1), y_predict) + 19*ydiff0 - 5*ydiff1 + ydiff2)*stepsize/24;YMat(i + 1, :) = y_corrector.';endendend。
题目:Matlab 解常微分方程的初值问题一、 设计目的:1、熟练掌握Matlab 的基本编程方法,及其编程风格。
2、熟练掌握Matlab 常用函数的使用。
3、与本专业相关知识相结合,掌握其在程序开发中的应用方法 以及和word 、C 语言等接口方法。
4、通过计算机数值求解的方式来加深微分方程解的理解。
5、熟悉初等方法可获得解析解之外的数值近似解的求解方法,提 高对差分格式的认识和离散化分析问题的技巧,加深对理论课程的学习和理解,为数学专业和信息与计算科学专业其他后继课程的学习打好基础。
二、 设计内容:已知一个三阶微分方程:012'''2'''=+⎪⎭⎫ ⎝⎛--y y y yy ,利用matlab软件求这个三阶微分方程在初值 ()()()100010'''-===y y y 下的解。
原三阶微分方程可化为:y y y y y '''2'''12-⎪⎭⎫ ⎝⎛-=令y y y y yy '''321=== 则原三阶微分方程可化为微分方程组213123213212'''y y y y y y y y y -⎪⎭⎫ ⎝⎛-=== 在初值 ()()()100010'''-===y y y 下的解。
三、 程序流程:四、程序代码:%编写函数文件rigid.mfunction 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.00010.00010.00020.00020.00050.0007………………0.93831.06651.19471.29181.38891.4860………………6.29166.29226.29286.29346.29406.29476.由于数据太多,这里只列举部分%绘制解的曲线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未加图形标注时的图012345678910Solution of Rigid Equationtime Ts o l u t i o n Y加了图形标注后的图输出结果[T,Y]中T为时间点组成的向量。
1用Matlab 求常微分方程(ODE)的初值问题(IVP)本节考虑一阶常微分方程000(,) ()u f t u t t Tu t u '=<≤⎧⎨=⎩ (1.1)的数值求解问题,包括算法公式及编程问题。
对一阶常微分方程组的初值问题010111120221220200120()(,,,,)(,,,,)() (,,,,)()m mm m m m m u t u u f t u u u u f t u u u u t u t t T u f t u u u u t u ⎧='=⎧⎪⎪'==⎪⎪<≤⎨⎨⎪⎪⎪⎪'==⎩⎩ (1.2)可以通过引入列向量0,,u u f化成类似(1.1)的形式000(,) ()u f t u t t Tu t u ⎧'=<≤⎪⎨=⎪⎩(1.3)其中1101122202120012()()(,,,,)()()(,,,,)(),,(,)()()(,,,,)m m m m mm u t u t f t u u u u t u t f t u u u u t u f t u u t u t f t u u u ⎛⎫⎛⎫⎛⎫⎪ ⎪ ⎪ ⎪ ⎪ ⎪=== ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭ (1.4)另外一个高阶常微分方程的初值问题()(1)0(1)(1)000000(,,,,,) (),(),,()m m m m u f t u u u u t t Tu t u u t u u t u ---'''⎫=<≤⎪⎬''===⎪⎭(1.5)也可以通过变换(1)123,,,,m m u u u u u u u u -'''==== 化成维微分方程组:1223112(,,,,)m mmm u u u u u uu f t u u u -'=⎧⎪'=⎪⎪⎨⎪'=⎪'=⎪⎩(1.6)我们在设计算法时通常先对一维一阶常微分方程(1.1)进行,然后再将这个算法写成适合求解(1.3)的向量形式,并以向量形式来进行编程。
matlab用经典龙格库塔法求微分方程组初值问题程序经典龙格-库塔法是一种数值求解常微分方程的方法。
以下是一个使用MATLAB实现经典龙格-库塔法求解微分方程组的示例代码:```matlabfunction [t, y] = runge_kutta(f, y0, tspan, N)% f: 微分方程右边的函数句柄% y0: 初始值% tspan: 时间范围 [t0, tf]% N: 步数t = linspace(tspan(1), tspan(2), N+1); % 时间向量y = zeros(size(t)); % 初始化解向量y(1) = y0; % 设置初始值for i = 1:N% 计算四个点上的值k1 = f(t(i), y(i));k2 = f(t(i) + h/2, y(i) + h/2k1);k3 = f(t(i) + h/2, y(i) + h/2k2);k4 = f(t(i) + h, y(i) + hk3);% 更新解向量y(i+1) = y(i) + h/6(k1 + 2k2 + 2k3 + k4);endend```在上述代码中,我们定义了一个名为 `runge_kutta` 的函数,该函数接受微分方程右边的函数句柄 `f`、初始值 `y0`、时间范围 `tspan` 和步数 `N` 作为输入,并返回时间向量 `t` 和解向量 `y`。
在函数内部,我们首先生成时间向量 `t`,然后初始化解向量 `y`,并设置初始值 `y0`。
接下来,我们使用一个循环来迭代计算每个时间点上的值,并使用龙格-库塔公式更新解向量。
最后,我们返回时间向量 `t` 和解向量 `y`。
要使用该函数求解微分方程组,可以按照以下步骤进行:1. 定义微分方程右边的函数句柄 `f`,该函数接受时间和解向量作为输入,并返回微分方程的右侧值。
2. 定义初始值 `y0`。
3. 定义时间范围 `tspan`,该向量包含时间范围的起始和终止值。
常微分方程的初值问题1.DEEuler 用欧拉法求一阶常微分方程的数值解function y = DEEuler(f, h,a,b,y0,varvec)format long;N = (b-a)/h;y = zeros(N+1,1);x = a:h:b;y(1) = y0;for i=2:N+1y(i) = y(i-1)+h*Funval(f,varvec,[x(i-1), y(i-1)]); endformat short;2.DEimpEuler 用隐式欧拉法求一阶常微分方程的数值解function y = DEimpEuler(f, h,a,b,y0,varvec)format long;N = (b-a)/h;y = zeros(N+1,1);y(1) = y0;x = a:h:b;var = findsym(f);for i=2:N+1fx = Funval(f,var(1),x(i));gx = y(i-1)+h*fx - varvec(2);y(i) = NewtonRoot(gx,-10,10,eps);endformat short;3.DEModifEuler 用改进欧拉法求一阶常微分方程的数值解function y = DEModifEuler(f, h,a,b,y0,varvec)format long;N = (b-a)/h;y = zeros(N+1,1);y(1) = y0;x = a:h:b;var = findsym(f);for i=2:N+1v1 = Funval(f,varvec,[x(i-1) y(i-1)]);t = y(i-1) + h*v1;v2 = Funval(f,varvec,[x(i) t]);y(i) = y(i-1)+h*(v1+v2)/2;endformat short;4.DELGKT2_mid 用中点法求一阶常微分方程的数值解function y = DELGKT2_mid(f, h,a,b,y0,varvec)format long;N = (b-a)/h;y = zeros(N+1,1);y(1) = y0;x = a:h:b;var = findsym(f);for i=2:N+1v1 = Funval(f,varvec,[x(i-1) y(i-1)]);t = y(i-1) + h*v1/2;v2 = Funval(f,varvec,[x(i)+h/2 t]);y(i) = y(i-1)+h*v2;endformat short;5.DELGKT2_suen 用休恩法求一阶常微分方程的数值解function y = DELGKT2_suen(f, h,a,b,y0,varvec)format long;N = uint16((b-a)/h);y = zeros(N+1,1);y(1) = y0;x = a:h:b;var = findsym(f);for i=2:N+1K1 = Funval(f,varvec,[x(i-1) y(i-1)]);K2 = Funval(f,varvec,[x(i-1)+2*h/3 y(i-1)+K1*2*h/3]); y(i) = y(i-1)+h*(K1+3*K2)/4;endformat short;6.DELGKT3_suen 用休恩三阶法求一阶常微分方程的数值解function y = DELGKT3_suen(f, h,a,b,y0,varvec)format long;N = (b-a)/h;y = zeros(N+1,1);y(1) = y0;x = a:h:b;var = findsym(f);for i=2:N+1K1 = Funval(f,varvec,[x(i-1) y(i-1)]);K2 = Funval(f,varvec,[x(i-1)+h/3 y(i-1)+K1*h/3]);K3 = Funval(f,varvec,[x(i-1)+2*h/3 y(i-1)+K2*2*h/3]); y(i) = y(i-1)+h*(K1+3*K3)/4;endformat short;7.DELGKT3_kuta 用库塔三阶法求一阶常微分方程的数值解function y = DELGKT3_kuta(f, h,a,b,y0,varvec)format long;N = (b-a)/h;y = zeros(N+1,1);y(1) = y0;x = a:h:b;var = findsym(f);for i=2:N+1K1 = Funval(f,varvec,[x(i-1) y(i-1)]);K2 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]);K3 = Funval(f,varvec,[x(i-1)+h y(i-1)-h*K1+K2*2*h]);y(i) = y(i-1)+h*(K1+4*K2+K3)/6;endformat short;8.DELGKT4_lungkuta 用经典龙格-库塔法求一阶常微分方程的数值解function y = DELGKT4_lungkuta(f, h,a,b,y0,varvec)format long;N = (b-a)/h;y = zeros(N+1,1);y(1) = y0;x = a:h:b;var = findsym(f);for i=2:N+1K1 = Funval(f,varvec,[x(i-1) y(i-1)]);K2 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]);K3 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K2*h/2]);K4 = Funval(f,varvec,[x(i-1)+h y(i-1)+h*K3]);y(i) = y(i-1)+h*(K1+2*K2+2*K3+K4)/6;endformat short;9.DELGKT4_jer 用基尔法求一阶常微分方程的数值解function y = DELGKT4_jer(f, h,a,b,y0,varvec)format long;N = (b-a)/h;C2 = sqrt(2);y = zeros(N+1,1);y(1) = y0;x = a:h:b;var = findsym(f);for i=2:N+1K1 = Funval(f,varvec,[x(i-1) y(i-1)]);K2 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]);K3 = Funval(f,varvec,[x(i-1)+h/2y(i-1)+K1*h*(C2-1)/2+K2*h*(2-C2)/2]);K4 = Funval(f,varvec,[x(i-1)+h y(i-1)-K2*h*C2/2+K3*h*(2+C2)/2]); y(i) = y(i-1)+h*(K1+(2-C2)*K2+(2+C2)*K3+K4)/6;endformat short;10.DELGKT4_qt 用变形龙格-库塔法求一阶常微分方程的数值解function y = DELGKT4_qt(f, h,a,b,y0,varvec)format long;N = (b-a)/h;y = zeros(N+1,1);y(1) = y0;x = a:h:b;var = findsym(f);for i=2:N+1K1 = Funval(f,varvec,[x(i-1) y(i-1)]);K2 = Funval(f,varvec,[x(i-1)+h/3 y(i-1)+K1*h/3]);K3 = Funval(f,varvec,[x(i-1)+2*h/3 y(i-1)-K1*h/3+K2*h]);K4 = Funval(f,varvec,[x(i-1)+h y(i-1)+h*(K1-K2+K3)]);y(i) = y(i-1)+h*(K1+3*K2+3*K3+K4)/8;endformat short;11.DELSBRK 用罗赛布诺克半隐式法求一阶常微分方程的数值解function y = DELSBRK(f, h,a,b,y0,varvec)format long;a1 = 1.40824829;a2 = 0.59175171;c1 = 0.17378667;c2 = c1;w1 = -0.41315432;w2 = 1.41315432;N = (b-a)/h;y = zeros(N+1,1);y(1) = y0;x = a:h:b;var = findsym(f);dy = diff(f, varvec(2));for i=2:N+1f1 = Funval(f,varvec,[x(i-1) y(i-1)]);dy1 = Funval(dy,varvec,[x(i-1) y(i-1)]);k1 = h*f1/(1-h*a1*dy1);dy2 = Funval(dy,varvec,[x(i-1)+c1*h y(i-1)+c2*k1]);f2 = Funval(f,varvec,[x(i-1)+c1*h y(i-1)+c2*k1]);k2 = h*f2/(1-h*a2*dy2);y(i) = y(i-1)+w1*k1+w2*k2;endformat short;12.DEMS 用默森单步法求一阶常微分方程的数值解function y = DEMS(f, h,a,b,y0,varvec)format long;N = (b-a)/h;y = zeros(N+1,1);y(1) = y0;x = a:h:b;var = findsym(f);for i=2:N+1fy = Funval(f,varvec,[x(i-1) y(i-1)]);z1 = y(i-1)+h*fy/3;fy1 = Funval(f,varvec,[x(i-1)+h/3 z1]);z2 = z1+h*(fy1-fy)/6;fy2 = Funval(f,varvec,[x(i-1)+h/3 z2]);z3 = z2+3*h*(fy2 - 4*fy1/9-fy/9)/8;fy3 = Funval(f,varvec ,[x(i-1)+h/2 z3]);z4 = z3+2*h*(fy3- 15*fy2/16 + 3*fy/16);fy4 = Funval(f,varvec ,[x(i-1)+h z4]);y(i) = z4+h*(fy4 - 8*fy3 + 9*fy2 - 2*fy)/6;endformat short;13.DEMiren 用米尔恩法求一阶常微分方程的数值解function y = DEMiren(f, h,a,b,y0,varvec)format long;N = (b-a)/h;y = zeros(N+1,1);x = a:h:b;y(1) = y0;y(2) = y0 + h*Funval(f,varvec, [x(1) y(1)]);y(3) = y(2) + h*Funval(f,varvec, [x(2) y(2)]);y(4) = y(3) + h*Funval(f,varvec, [x(3) y(3)]);var = findsym(f);for i=5:N+1y(i) = y(i-4)+4*h*(2*Funval(f,varvec,[x(i-1), y(i-1)]) - ... Funval(f,varvec,[x(i-2), y(i-2)]) + ...2*Funval(f,varvec,[x(i-3), y(i-3)]))/3; endformat short;14.DEYDS 用亚当斯法求一阶常微分方程的数值解function y = DEYDS(f, h,a,b,y0,varvec)format long;N = (b-a)/h;y = zeros(N+1,1);x = a:h:b;y(1) = y0;y(2) = y0 + h*Funval(f,varvec, [x(1) y(1)]);y(3) = y(2) + h*Funval(f,varvec, [x(2) y(2)]);y(4) = y(3) + h*Funval(f,varvec, [x(3) y(3)]);y(5) = y(4) + h*Funval(f,varvec, [x(4) y(4)]);var = findsym(f);for i=6:N+1y(i) = y(i-1)+h*(55*Funval(f,varvec,[x(i-1), y(i-1)])/24 - ... 59*Funval(f,varvec,[x(i-2), y(i-2)])/24 + ... 37*Funval(f,varvec,[x(i-3), y(i-3)])/24 - ... 9*Funval(f,varvec,[x(i-4), y(i-4)])/24);endformat short;15.DEYCJZ_mid 用中点-梯形预测校正法求一阶常微分方程的数值解function y = DEYCJZ_mid(f, h,a,b,y0,varvec,type,s)format long;N = (b-a)/h;y = zeros(N+1,1);x = a:h:b;y(1) = y0;y(2) = y0+h*Funval(f,varvec,[x(1) y(1)]);var = findsym(f);if type == 1for i=3:N+1v1 = Funval(f,varvec,[x(i-1) y(i-1)]);t = y(i-2) + 2*h*v1;v2 = Funval(f,varvec,[x(i) t]);y(i) = y(i-1)+h*(v1+v2)/2;endelseif type == 2for i=3:N+1v1 = Funval(f,varvec,[x(i-1) y(i-1)]);t = y(i-2) + 2*h*v1;v2 = Funval(f,varvec,[x(i) t]);for l=1:sy(i) = y(i-1)+h*(v2 + v1)/2;v2 = Funval(f,varvec,[x(i) y(i)]);endendelsep0 = 0;c = 0;fnl = Funval(f,varvec,[x(1) y(1)]);for i = 3:N+1v1 = Funval(f,varvec,[x(i-1) y(i-1)]);p = y(i-2)+2*h*v1;M = p - 0.8*(p0 - c);F = Funval(f , varvec, [x(i) ,M]);fn = y(i-1) + h*(fnl + F)/2;y(i) = fn - 0.2*( p - fn);c = fn;fnl = fn;p0 = p;endendendformat short;16.DEYCJZ_adms 用阿达姆斯预测校正法求一阶常微分方程的数值解function y = DEYCJZ_adms(f, h,a,b,y0,varvec,type) format long;N = (b-a)/h;y = zeros(N+1,1);x = a:h:b;y(1) = y0;y(2) = y0+h*Funval(f,varvec,[x(1) y(1)]);for i=3:N+1v1 = Funval(f,varvec,[x(i-2) y(i-2)]);v2 = Funval(f,varvec,[x(i-1) y(i-1)]);t = y(i-1) + h*(3*v2-v1)/2;v3 = Funval(f,varvec,[x(i) t]);y(i) = y(i-1)+h*(v2+v3)/2;endformat short;17.DEYCJZ_adms2 用密伦预测校正法求一阶常微分方程的数值解function y = DEYCJZ_adms2(f, h,a,b,y0,varvec,type) format long;N = (b-a)/h;y = zeros(N+1,1);x = a:h:b;y(1) = y0;y(2) = y0+h*Funval(f,varvec,[x(1) y(1)]);y(3) = y(2)+h*Funval(f,varvec,[x(2) y(2)]);y(4) = y(3)+h*Funval(f,varvec,[x(3) y(3)]);for i=5:N+1v2 = Funval(f,varvec,[x(i-3) y(i-3)]);v3 = Funval(f,varvec,[x(i-2) y(i-2)]);v4 = Funval(f,varvec,[x(i-1) y(i-1)]);t = y(i-1) + h*(55*v4 - 59*v3 + 37*v2 - 9*v1)/24;ft = Funval(f,varvec,[x(i) t]);y(i) = y(i-1)+h*(9*ft+19*v4-5*v3+v2)/24;endformat short;18.DEYCJZ_ yds 用亚当斯预测校正法求一阶常微分方程的数值解function y = DEYCJZ_yds (f, h,a,b,y0,varvec)format long;N = (b-a)/h;y = zeros(N+1,1);x = a:h:b;y(1) = y0;y(2) = y0+h*Funval(f,varvec,[x(1) y(1)]);y(3) = y(2)+h*Funval(f,varvec,[x(2) y(2)]);y(4) = y(3)+h*Funval(f,varvec,[x(3) y(3)]);for i=5:N+1v1 = Funval(f,varvec,[x(i-4) y(i-4)]);v2 = Funval(f,varvec,[x(i-3) y(i-3)]);v3 = Funval(f,varvec,[x(i-2) y(i-2)]);v4 = Funval(f,varvec,[x(i-1) y(i-1)]);t = y(i-1) + h*(55*v4 - 59*v3 + 37*v2 - 9*v1)/24;ft = Funval(f,varvec,[x(i) t]);y(i) = y(i-1)+h*(9*ft+19*v4-5*v3+v2)/24;endformat short;19.DEYCJZ_ myds 用修正的亚当斯预测校正法求一阶常微分方程的数值解function y = DEYCJZ_myds(f, h,a,b,y0,varvec,type)format long;N = (b-a)/h;y = zeros(N+1,1);x = a:h:b;y(1) = y0;y(2) = y0+h*Funval(f,varvec,[x(1) y(1)]);y(3) = y(2)+h*Funval(f,varvec,[x(2) y(2)]);y(4) = y(3)+h*Funval(f,varvec,[x(3) y(3)]);p0 = 0;c = 0;for i=5:N+1v2 = Funval(f,varvec,[x(i-3) y(i-3)]);v3 = Funval(f,varvec,[x(i-2) y(i-2)]);v4 = Funval(f,varvec,[x(i-1) y(i-1)]);p = y(i-1) + h*(55*v4 - 59*v3 + 37*v2 - 9*v1)/24;M = p + 251*(c - p0)/270;F = Funval(f,varvec,[x(i) M]);c = y(i-1) + h*(9*F + 19*v4 - 5*v3 + v2)/24;y(i) = c+h*19*(p - c)/270;p0 = p;endformat short;20.DEYCJZ_hm 用汉明预测校正法求一阶常微分方程的数值解function y = DEYCJZ_hm(f, h,a,b,y0,varvec,type)format long;N = (b-a)/h;y = zeros(N+1,1);x = a:h:b;y(1) = y0;y(2) = y0+h*Funval(f,varvec,[x(1) y(1)]);y(3) = y(2)+h*Funval(f,varvec,[x(2) y(2)]);y(4) = y(3)+h*Funval(f,varvec,[x(3) y(3)]);if type == 1for i=5:N+1v1 = Funval(f,varvec,[x(i-3) y(i-3)]);v2 = Funval(f,varvec,[x(i-2) y(i-2)]);v3 = Funval(f,varvec,[x(i-1) y(i-1)]); t = y(i-4) + 4*h*(2*v3 - v2 + 2*v1)/3;ft = Funval(f,varvec,[x(i) t]);y(i) = (9*y(i-1) -y(i-3) +3*h*(2*v3 + ft-v2))/8;endelsep0 = 0;c = 0;for i = 5:N+1v1 = Funval(f,varvec,[x(i-3) y(i-3)]);v2 = Funval(f,varvec,[x(i-2) y(i-2)]);v3 = Funval(f,varvec,[x(i-1) y(i-1)]);p = y(i-4)+4*h*(2*v3 - v2 + 2*v1)/3;M = p - 112*(p0 - c)/121;F = Funval(f , varvec, [x(i) ,M]);c = (9*y(i-1) -y(i-3) +3*h*(2*v3 + F-v2))/8;y(i) = c + 9*( p - c)/121;p0 = p;endendformat short;21.DEWT 用外推法求一阶常微分方程的数值解function yy = DEWT(f,h,a,b,gama,y0,order,varvec)format long;ArrayH = [1;2;4;6;8;12;16;24;32;48;64;96];N = (b-a)/h;yy = zeros(N+1,1);for i = 2:N+1dh = h;s = zeros(order,1);for j=1:orderdh = h/ArrayH(j);tmpY = DELGKT2_suen(f,dh,a,a+(i-1)*h,y0,varvec);s(j) = tmpY((i-1)*ArrayH(j)+1);endtmpS = zeros(order,1);for j=1:order-1for k=(j+1):ordertmpS(k) = s(k)+(s(k)-s(k-1))/((ArrayH(k)/ArrayH(j))^gama-1); ends(1:(order-j)) = tmpS((j+1):order);endyy(i) = tmpS(order);endformat short;22.DEWT_glg 用格拉格外推法求一阶常微分方程的数值解function yy = DEWT_glg(f,h,a,b,gama,y0,order,varvec)%一阶常微分方程的一般表达式的右端函数:f%积分步长:h%自变量取值下限:a%自变量取值上限:b%外推参数,参考外推公式:gama%函数初值:y0%外推阶数:order%常微分方程的变量组:varvecArrayH = [1;2;4;6;8;12;16;24;32;48;64;96];N = (b-a)/h;yy = zeros(N+1,1);for i = 2:N+1dh = h;s = zeros(order,1);for j=1:orderdh = h/ArrayH(j);NI = (i-1)* ArrayH(j);tmpY = zeros(NI+2,1);tmpY(1) = y0;tmpY(2) = y0 + dh*Funval(f,varvec,[a+dh y0]);for m=3:NI+2tmpY(m) = tmpY(m-2)+ 2*dh*Funval(f,varvec,[a+(m-1)*dh tmpY(m-1)]);ends(j) = tmpY(NI)/4+ tmpY(NI+1)/2+ tmpY(NI+2)/4;endtmpS = zeros(order,1);for j=1:order-1for k=(j+1):ordertmpS(k) = s(k)+(s(k)-s(k-1))/((ArrayH(k)/ArrayH(j))^gama-1); ends(1:(order-j)) = tmpS((j+1):order);endyy(i) = tmpS(order);end。