数学实验“欧拉数值算法,Runge-Kutta数值算法”实验报告(内含matlab程序)
数学实验四
实验四1415100076 数学141 陈田实验项目名称：常微分方程数值解实验时间：2016-4-29、2016-5-6实验地点：理学楼525实验目的：1、用Matlab软件掌握求微分方程数值解的欧拉方法和龙格-库塔方法；2、通过实例学习用微分方程模型解决简化的实际问题。
实验内容：1、编写用向前欧拉公式和改进欧拉公式求微分方程数值解的程序。
function [y1,y2] = heuler( x0,y0,b,h)%x0,y0:初始条件%x0,b:求值区间%h:步长x=x0:h:b;n=(b-x0)/h+1;y1(1)=y0;y2(1)=y0;for i=2:ny1(i)=y1(i-1)+h.*ceuler(x(i-1),y1(i-1)); %向前欧拉公式endfor j=2:ny2(j)=y2(j-1)+h/2.*(ceuler(x(j-1),y2(j-1))+ceuler(x(j),y2(j-1)+h.*ceu ler(x(j-1),y2(j-1)))); %改进欧拉公式endplot(x,y1,'g+',x,y2,'r+',x,3.*exp(x)-2.*x-2,'b')end2、教材120页第1题a；(提示：原微分方程应转化为参数方程，以适应MATLAB的龙格-库塔函数调用格式。
可令参数t=x，则dx/dt=1,dy/dt=dy/dx=y+2t，然后建立参数方程的函数M文件。
) function d = canshu( t,x ) %建立参数方程的函数M文件d=[1;x(2)+2*t];endM文件[y1,y2] = heuler( 0,1,1,0.1);t0=0:0.1:1;x0=[0,1];[t,x]=ode45(@canshu,t0,x0);[t,x]plot(x(:,1),x(:,2),'kd',t0,y1,'g+',t0,y2,'r+',t0,3.*exp(t0)-2.*t0-2,' b')运行结果>> rkfbans =0 0 1.00000.1000 0.1000 1.11550.2000 0.2000 1.26420.3000 0.3000 1.44960.4000 0.4000 1.67550.5000 0.5000 1.94620.6000 0.6000 2.26640.7000 0.7000 2.64130.8000 0.8000 3.07660.9000 0.9000 3.57881.0000 1.0000 4.15483、教材122页第5题，不做题目要求中的求解析解及数值解与解析解作比较。
《数值分析》第五章实验报告
1.900 11.7479965 2.000 15.3982357 则有 i 1 5 6 9 10 ti 1.1 1.5 1.6 1.9 2.0 wi 0.2718282 3.1874451 4.6208178 11.7479965 15.3982357 y(ti) 0.345920 3.96767 5.70296 14.3231 18.6831
b)c)d)类似进行即可
EXERCISE SET 5.9 P322 2、方程组的 Runge-Kutta 算法 a) y' '2 y' y te t ,0 t 1, y(0) y' (0) 0, h 0.1
t
设 u1 (t ) y(t ), u2 (t ) y (t ) ，则将方程转换为方程组
'
-5-
u1' (t ) u2 (t )
' u2 (t ) 2u2 (t ) u1 (t ) t (et 1)
初始条件为
u1 (0) 0, u2 (0) 0
编写 MATLAB 程序 function[t,y] = Runge_Kutta4s(ydot_fun,t0,y0,h,N) %标准四阶Runge_Kutta公式，其中， %ydot_fun为一阶微分方程的函数； %t0为初始点； %y0为初始向量（列向量） ； %h为区间步长； %N为区间的个数； %t为Tn构成的向量； %y为Yn构成的矩阵。 t = zeros(1,N+1);y = zeros(length(y0),N+1); t(1) = t0;y(:,1) = y0; for n = 1 :N t(n+1) = t(n) + h; k1 = h * feval(ydot_fun,t(n),y(:,n)); k2 = h * feval(ydot_fun,t(n)+1/2 * h,y(:,n)+1/2 * k1); k3 = h * feval(ydot_fun,t(n)+1/2 * h,y(:,n)+1/2 * k2); k4 = h * feval(ydot_fun,t(n)+h,y(:,n)+k3); y(:,n+1) = y(:,n) + 1/6 * (k1 + k2 + k3 + k4); end 运行后有 >> odefun = inline('[y(2);2*y(2)-y(1)+t*(exp(t)-1)]','t','y'); >> [t,y] = Runge_Kutta4s(odefun,0,[0;0],0.1,10) t= Columns 1 through 9 0 0.8000 Columns 10 through 11 0.9000 1.0000 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000
Runge-Kutta报告
实验报告——四级四阶Runge-Kutta 公式1110500158 刘雪婷一 实验问题从上一个报告中知Euler 法和梯形法是最简单的求解常微分方程离散点值的方法。
但这两种方法的精确度都比较低。
在解微分方程的时候为了使精度更高，我们使用了Taylor 级数法，通过增加展开项数来提高方法的阶。
但是，这种方法计算量大，所以我们用一种间接使用Taylor 展开式来获取高阶的方法。
这种方法就是Runge-Kutta 法，既可以减少计算量，又可以提高精度。
二 实验方法{ du dt =ℱ(t,u)u (t 0)=u 0* 对*式左右两边积分，得到： u (t n+1)−u (t n )=∫ℱ(t,u(t))dt t n+1t n 积分中值定理1 ⇒ h ℱ(t n +θh,u(t n +θh))我们用［t n ,t n+1］中若干个点的线性组合来近似θ。
u n+1=u (t n )+h ∑ωi k i s i=1#k i =ℱ(t n +hαi ,u (t n )+h ∑βi,j k j i−1j=1) ,i ≥1 k 1=ℱ(t n ,u(t n )) 将#式两边在t n 处Taylor 展开，然后按照h 的幂排列，让同幂的两项系数相等，以此来确定ωi αj βij 。
利用上述方法，我们得到了古典Rung-Kutta 公式（四级四阶）： y n+1=y n +h 6(K 1+2K 2+2K 3+K 4)K 1=ℱ(x n ,y n )K 2=ℱ(x n +h 2,y n +h 2K 1 K 3=ℱ(x n +h 2,y n +h 2K 2 K 4=ℱ(x n +h,y n +hK 3 三 初始量{dy dx =x 2−y y (0)=10≤x ≤1 h=0.1四 结论可以看到，RK 方法得到的结果精确到小数点后5位。
五 结论分析当h=0.2时，结果精确到小数点后4位。
当h=0.01时，结果精确到了小数点后9位。
数值分析第七章实验报告7
贵州师范大学数学与计算机科学学院学生实验报告课程名称： 数值分析 班级： 实验日期： 2011年12月14日 学号： 姓名： 指导教师： 实验成绩： 一、实验名称实验六: 常微分方程初值问题数值解法 二、实验目的及要求1. 让学生掌握用Euler 法, Runge-Kutta 法求解常微分方程初值问题.2. 培养Matlab 编程与上机调试能力.三、实验环境每人一台计算机，要求安装Windows XP 操作系统，Microsoft office2003、MATLAB6.5(或7.0). 四、实验内容1. 取步长h=0.1,0.05,0.01,,用Euler 法及经典4阶Runge-Kutta 法求解初值问题⎩⎨⎧=≤≤++-=1)0()10(2222'y t t t y y 要求:1) 画出准确解(准确解22t e y t +=-)的曲线,近似解折线;2) 把节点0.1和0.5上的精确解与近似解比较,观察误差变化情况.2. 用 Euler 法,隐式Euler 法和经典4阶R-K 法取不同步长解初值问题⎪⎩⎪⎨⎧=∈-=21)0(],1,0[,50'y x y y 并画出曲线观察稳定性. 注：题1必须写实验报告五、算法描述及实验步骤 Euler 法：输入 000),(,,,),,(y a x x h b a y x f = 输出 Euler 解y 步1 ),,2,1(;m n h n a x ha b m n =⨯+=-⇐步2 对1,,2,1,0-=m n 执行),(1n n n n y x f h y y ⨯+⇐+步3 输出T m y y y y ),,,(21 = 经典4阶R-K 法：输入 000),(,,,),,(y a x x h b a y x f = 输出 4阶R-K 解y 步1 ),,2,1(;m n h n a x ha b m n =⨯+=-⇐步2 对1,,2,1,0-=m n 执行),(1n n y x f K ⇐，)5.0,(15.02hK y x f K n n +⇐+， )5.0,(25.03hK y x f K n n +⇐+，),(314hK y x f K n n +⇐+)22(643211K K K K h y y n n ++++⇐+步3 输出T m y y y y ),,,(21 = 六、调试过程及实验结果 分别取h=0.1,h=0.01 步长h=0.1>> a=0;b=1;h=0.1;y0=1; >> Y1=RK('fun',a,b,y0,h); >> y1=[1,Y1];>> Y2=Euler('fun',a,b,y0,h); >> y2=[1,Y2]; >> t=[0:0.1:1]; >> plot(t,y1,'dk-') >> hold on>> plot(t,y2,'.r--') >> t=0:0.01:1; >> y=exp(-2.*t)+t.^2; >> hold on >> plot(t,y,':')>> legend('y1','y2','y',2)>> title('准确解曲线与近似解折线')>> text(0.15,0.775,'\leftarrowy1=R-K 法') >> text(0.7,0.67,'\leftarrowy2=Euler 法') >> text(0.48,0.78,'y=准确解曲线\rightarrow')>> y=exp(-2.*t)+t.^2; >> t=0.1;>> y;>> Y1(1);>> Y2(2);>> e1=y-Y1(1)>> y=exp(-2.*t)+t.^2; >> t=0.5;>> y;>> Y1(5);>> Y2(5);>> e1=y-Y1(5)e1 =-1.1609e-005>> e2=y-Y2(5)e2 =0.0738e1 =-4.2469e-006>> e2=y-Y2(1)e2 =0.0287步长h=0.01>> a=0;b=1;h=0.01;y0=1;>> Y1=RK('fun',a,b,y0,h); >> y1=[1,Y1];>> Y2=Euler('fun',a,b,y0,h); >> y2=[1,Y2];>> y=exp(-2.*t)+t.^2;>> t=0:0.01:1;>> plot(t,y1,'k-')>> hold on>> plot(t,y2,'r--')>> hold on>> plot(t,y,':')>> t=0.1;>> y=exp(-2.*t)+t.^2;>> y;>> Y1(10);>> Y2(10);>> e1=y-Y1(10)e1 =-3.7457e-010>> e2=y-Y2(10)e2 =0.0026>> t=0.5;>> y=exp(-2.*t)+t.^2;>> y;>> Y1(50);>> Y2(50);>> e1=y-Y1(50)e1 =-1.0308e-009>> e2=y-Y2(50)e2 =0.0069七、总结用1阶的Euler法解初值问题时，若步长过大的话，误差将会较大，其解不可靠，只有控制步长尽量小，在一定误差范围内才可用，大的话精度不高。
matlab-欧拉方法和龙格库塔方法的小实例
题一：a)y’=y+2x ， 欧拉方法：112()2n n h y y k k +=++，12n n k y x =+，2112()n n k y hk x +=++； 龙格-库塔方法：11234(22)6n n h y y k k k k +=++++，12n n k y x =+，12222n n k h k y h x ⎛⎫=+++ ⎪⎝⎭，23222n n k h k y h x ⎛⎫=+++ ⎪⎝⎭，432()n n k y hk x h =+++ 精确解：y=3e x -2x-2。
以步长h=0.1 在0<=x<=1内的计算结果如下所示：0.1000 1.1150 1.1155 1.11550.2000 1.2631 1.2642 1.26420.3000 1.4477 1.4496 1.44960.4000 1.6727 1.6755 1.67550.5000 1.9423 1.9462 1.94620.6000 2.2613 2.2664 2.26640.7000 2.6347 2.6413 2.64130.8000 3.0684 3.0766 3.07660.9000 3.5685 3.5788 3.57881.0000 4.1422 4.1548 4.1548b)文案 编辑词条B 添加义项 ?文案，原指放书的桌子，后来指在桌子上写字的人。
现在指的是公司或企业中从事文字工作的职位，就是以文字来表现已经制定的创意策略。
文案它不同于设计师用画面或其他手段的表现手法，它是一个与广告创意先后相继的表现的过程、发展的过程、深化的过程，多存在于广告公司，企业宣传，新闻策划等。
基本信息中文名称文案外文名称Copy目录1发展历程2主要工作3分类构成4基本要求5工作范围6文案写法7实际应用折叠编辑本段发展历程汉字"文案"(wén àn)是指古代官衙中掌管档案、负责起草文书的幕友,亦指官署中的公文、书信等;在现代，文案的称呼主要用在商业领域，其意义与中国古代所说的文案是有区别的。
数值分析Runge现象计算实验
数值分析实验报告（02）一、实验目的通过上机绘制Runge 函数图像，理解高次插值的病态性质。
二、实验内容在区间[-1,1]上分别取n=10,n=20用两组等距节点对龙格(Runge)函数21()125f x x =+作多项式插值，对每个n 值分别画出()f x 和插值函数的图形。
三、编程思路（相关背景知识、算法步骤、流程图、伪代码）四、程序代码（Matlab 或C 语言的程序代码）function yt=Untitled8(x,y,xt)%UNTITLED5 ´Ë´¦ÏÔʾÓйش˺¯ÊýµÄÕªÒª% ´Ë´¦ÏÔʾÏêϸ˵Ã÷n=length(x);ny=length(y);if n~=nyerror('²åÖµ½ÚµãxÓ뺯ÊýÖµy²»Ò»ÖÂ');endm=length(xt);yt=zeros(1,m);for k=1:nlk=ones(1,m);for j=1:nif j~=klk=lk.*(xt-x(j))/(x(k)-x(j));endend ;yt=yt+y(k)*lk;endn=input('n=');x=linspace(-1,1,n);y=1./(1+25.*x.^2);xf=linspace(-1,1,100);yf=1./(1+25.*xf.^2)xl=xf;yl=Untitled8(x,y,xf);plot(xf,yf,'-b',xl,yl,'-r')五、数值结果及分析（数值运行结果及对结果的分析）当n=10时当n=20六、实验体会（计算中出现的问题，解决方法，实验体会）出现符号错误，代码函数变量不明重新输入，查询错误，找到并改正编码需要认真仔细，一定要头脑清晰，避免出现一些低级错误。
龙格库塔实验报告
一、实验背景常微分方程（ODE）在自然科学、工程技术等领域中具有广泛的应用。
然而，许多微分方程无法得到精确解析解，因此需要借助数值方法进行求解。
龙格-库塔（Runge-Kutta）方法是一种常用的数值求解常微分方程的方法，具有精度高、稳定性好等优点。
本实验旨在通过编写程序，实现四阶龙格-库塔方法，并验证其在求解常微分方程中的有效性和准确性。
二、实验目的1. 理解四阶龙格-库塔方法的基本原理和计算步骤。
2. 编写程序实现四阶龙格-库塔方法。
3. 选取典型常微分方程，验证四阶龙格-库塔方法的求解精度和稳定性。
三、实验原理四阶龙格-库塔方法是一种基于泰勒级数展开的数值方法，其基本思想是将微分方程的解在某个区间内进行近似，并通过迭代计算得到近似解。
具体步骤如下：1. 初始化：给定初始条件y0，步长h，求解区间[a, b]。
2. 迭代计算：对于k=1, 2, ..., n（n为迭代次数），- 计算k1 = f(xk-1, yk-1)（f为微分方程的右端函数）；- 计算k2 = f(xk-1 + h/2, yk-1 + h/2 k1)；- 计算k3 = f(xk-1 + h/2, yk-1 + h/2 k2)；- 计算k4 = f(xk-1 + h, yk-1 + h k3)；- 更新y值：yk = yk-1 + (h/6) (k1 + 2k2 + 2k3 + k4)；- 更新x值：xk = xk-1 + h；3. 输出结果：输出最终的近似解y(n)。
四、实验步骤1. 编写程序实现四阶龙格-库塔方法。
2. 选取典型常微分方程，如：- y' = -y，初始条件y(0) = 1，求解区间[0, 2π]；- y' = y^2，初始条件y(0) = 1，求解区间[0, 1]。
3. 对每个常微分方程，设置不同的步长h和迭代次数n，分别计算近似解y(n)。
4. 将计算得到的近似解与解析解进行比较，分析四阶龙格-库塔方法的精度和稳定性。
Runge-Kutta法求微分方程数值解及其Matlab实现开题报告范文
毕业论文（设计）材料题目：Runge-Kutta法求微分方程数值解及其Matlab实现学生姓名：程晓曦学生学号：1005020103系别：数学与计算科学系专业：数学与应用数学届别：2010届指导教师：李宁填写说明1、本材料包括淮南师范学院本科毕业论文（设计）任务书、开题报告以及毕业论文（设计）评审表三部分内容。
2、本材料填写顺序依次为：（1）指导教师下达毕业论文（设计）任务书；（2）学生根据毕业论文（设计）任务书的要求，在文献查阅的基础上撰写开题报告，送交指导教师审阅并签字认可；（3）毕业论文（设计）工作后期，学生填写毕业论文（设计）主要内容，连同毕业论文（设计）全文一并送交指导教师审阅，指导教师根据学生实际完成的论文（设计）质量进行评价；（4）指导教师将此表连同学生毕业论文（设计）全文一并送交评阅教师评阅。
3、指导教师、评阅教师对学生毕业论文（设计）的成绩评定均采用百分制。
4、毕业论文（设计）答辩记录不包括在此表中。
一、毕业论文（设计）任务书二、毕业论文（设计）开题报告三、毕业论文（设计）评审表原文已完。
下文为附加文档，如不需要，下载后可以编辑删除，谢谢！施工组织设计本施工组织设计是本着“一流的质量、一流的工期、科学管理”来进行编制的。
编制时，我公司技术发展部、质检科以及项目部经过精心研究、合理组织、充分利用先进工艺，特制定本施工组织设计。
一、工程概况：西夏建材城生活区27#、30#住宅楼位于银川市新市区，橡胶厂对面。
本工程由宁夏燕宝房地产开发有限公司开发，银川市规划建筑设计院设计。
本工程耐火等级二级，屋面防水等级三级，地震防烈度为8度，设计使用年限50年。
本工程建筑面积:27#楼3824.75m2;30#楼3824.75 m2。
室内地坪±0.00以绝对标高1110.5 m为准，总长27#楼47.28m；30#楼47.28 m。
总宽27#楼14.26m；30#楼14.26 m。
设计室外地坪至檐口高度18.6 00m，呈长方形布置，东西向，三个单元。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
西京学院数学软件实验任务书实验二十四实验报告一、实验名称：欧拉数值算法（显式，隐式，欧拉预估-校正法），Runge-Kutta 数值算法。
二、实验目的：进一步熟悉欧拉数值算法（显式，隐式，欧拉预估-校正法），Runge-Kutta 数值算法。
三、实验要求：运用Matlab/C/C++/Java/Maple/Mathematica 等其中一种语言完成程序设计。
四、实验原理：1．欧拉数值算法（显式）：微分方程里最简单的方程形式莫过于一阶常微分方程的初值问题，即：0(,)()dy f x y a x b dx y a y ⎧=≤≤⎪⎨⎪=⎩其中a ，b 为常数。
因为其简单但又是求解其他方程的基础，所以发展了许多典型的解法。
所有算法中的f 就是代表上式中(,)f x y ，而y f 表示(,)f x y y ∂∂，x f 表示(,)f x y x∂∂。
简单欧拉法是一种单步递推算法。
简单欧拉法的公式如下所示：1(,)n n n n y y hf x y +=+简单欧拉法的算法过程介绍如下:给出自变量x 的定义域[,]a b ，初始值0y 及步长h 。
对0,1,()/k b a h =-，计算1(,)k k k k y y hf x y +=+2．欧拉数值算法（隐式）：隐式欧拉法也叫退欧拉法，隐式欧拉法的公式如下所示：111(,)n n n n y y hf x y +++=+隐式欧拉法是一阶精度的方法，比它精度高的公式是：111[(,)(,)]2n n n n n n h y y f x y f x y +++=++ 隐式欧拉的算法过程介绍如下。
给出自变量x 的定义域[,]a b ，初始值0y 及步长h 。
对0,1,()/k b a h =-，用牛顿法或其他方法求解方111(,)k k k k y y hf x y +++=+得出1k y +。
3．欧拉预估-校正法：改进欧拉法是一种二阶显式求解法，其计算公式如下所示：1[,(,)]22n n n n n n h h y y hf x y f x y +=+++11(,)[(,)(,)]2n n n n n n n n t y hf x y h y y f x y f x t ++=+⎧⎪⎨=++⎪⎩ 4．Runge-Kutta 数值算法：二阶龙格-库塔法有多种形式，除了改进的欧拉法外，还有中点法。
中点法计算公式为：1[,(,)]22n n n n n n h h y y hf x y f x y +=+++ 五、实验内容：%欧拉数值算法（显式）function [h,k,X,Y,P,REn]=Qeuler1(funfcn,x0,y0,b,n,tol) x=x0; h=(b-x)/n; X=zeros(n,1); y=y0;Y=zeros(n,1); k=1; X(k)=x; Y(k)=y';for k=2:n+1fxy=feval(funfcn,x,y);delta=norm(h*fxy,'inf');wucha=tol*max(norm(y,'inf'),1.0);if delta>=wuchax=x+h; y=y+h*fxy; X(k)=x;Y(k)=y';endplot(X,Y,'rp')gridlabel('自变量 X'), ylabel('因变量 Y')title('用向前欧拉（Euler ）公式计算dy/dx=f(x,y)，y(x0)=y0在[x0,b]上的数值解')endP=[X,Y];syms dy2,REn=0.5*dy2*h^2;%欧拉数值算法（隐式）function [X,Y,n,P]=Heuler1(funfcn,x0,b,y0,h,tol)n=fix((b-x0)/h); X=zeros(n+1,1); Y=zeros(n+1,1);k=1; X(k)=x0; Y(k,:)=y0; Y1(k,:)=y0;clc,x0,h,y0for i=2:n+1X(i)=x0+h; Y(i,:)=y0+h*feval(funfcn,x0,y0);Y1(i,:)=y0+h*feval(funfcn,X(i),Y(i,:));Wu=abs(Y1(i,:)-Y(i,:));while Wu>tolp=Y1(i,:);Y1(i,:)=y0+h*feval(funfcn,X(i),p);Y(i,:)=p;endx0=x0+h; y0=Y1(i,:);Y(i,:)=y0; plot(X,Y,'ro')grid onxlabel('自变量 X'), ylabel('因变量 Y')title('用向后欧拉公式计算dy/dx=f(x,y)，y(x0)=y0在[x0,b]上的数值解')endX=X(1:n+1); Y=Y(1:n+1,:); n=1:n+1; P=[n',X,Y]%欧拉预估-校正法function [H,X,Y,k,h,P]=gaiEuler(funfcn,x0,b,y0,tol) pow=1/3;if nargin < 5 | isempty(tol),tol = 1.e-6;end;if nargin < 6 | isempty(trace),trace = 0;end;x=x0; h=0.0078125*(b-x);y=y0(:);p=128; n=fix((b-x0)/h);H=zeros(p,1); X=zeros(p,1);Y=zeros(p,length(y)); k=1;X(k)=x; Y(k,:)=y';while (x<b)&(x+h>x)if x+h>bh=b-x;endfxy=feval(funfcn,x,y); fxy=fxy(:);delta=norm(h*fxy,'inf');wucha=tol*max(norm(y,'inf'),1.0);if delta<=wuchax=x+h; y1=y+h*fxy; fxy1=feval(funfcn,x,y1);fxy=fxy(:);y2=(fxy+fxy1)/2; y=y+h*y2; k=k+1;if k>length(X)X=[X;zeros(p,1)]; Y=[Y;zeros(p,length(y))];H=[H;zeros(p,1)];endH(k)=h;X(k)=x;Y(k,:)=y'; plot(X,Y,'mh'), gridxlabel('自变量 X'), ylabel('因变量 Y')title('用改进的欧拉公式计算dy/dx=f(x,y)，y(x0)=y0在[x0,b]上的数值解')endif delta~=0.0h=min(h*8,0.9*h*(wucha/delta)^pow);endendif (x<b)disp('Singularity likely.'), xendH=H(1:k); X=X(1:k); Y=Y(1:k,:); n=1:k; P=[n',H,X,Y] %Runge-Kutta数值算法function[k,X,Y,fxy,wch,wucha,P]=RK3(funfcn,fun,x0,b,C,y0,h)x=x0; y=y0;p=128;n=fix((b-x0)/h);fxy=zeros(p,1);wucha=zeros(p,1);wch=zeros(p,1);X=zeros(p,1);Y=zeros(p,length(y));k=1; X(k)=x; Y(k,:)=y';clc,x,h,yfor k=2:n+1x=x+h;a2=C(4);a3=C(5);b21=C(6);b31=C(7);b32=C(8);c1=C(1);c2=C(2);c3=C(3);x1=x+a2*h;x2=x+a3*h;k1=feval(funfcn,x,y);y1=y+b21*h*k1;k2=feval(funfcn,x1,y1);y2=y+b31*h*k1+b32*h*k2;k3=feval(funfcn,x2,y2);fxy(k)=feval(fun,x);y=y+h*(c1*k1+c2*k2+c3*k3);X(k)=x; Y(k,:)=y;k=k+1;plot(X,Y,'rp',X,fxy,'bo'),grid,xlabel('自变量 X'), ylabel('因变量 Y')legend('用三阶龙格-库塔方法计算dy/dx=f(x,y)，y(x0)=y0在[x0,b]上的数值解','y/dx=f(x,y)，y(x0)=y0的精确解y=f(x)') endfor k=2:n+1wucha(k)=norm(Y(k-1)-Y(k));wch(k)=norm(fxy(k)-Y(k));endX=X(1:k);Y=Y(1:k,:);fxy=fxy(1:k,:);n=1:k;wucha=wucha(1:k,:);wch=wch(1:k,:);P=[n',X,Y,fxy,wch,wucha];。