解微分方程欧拉法-R-K法及其MATLAB实例
龙哥库塔法or欧拉法求解微分方程matlab实现
龙哥库塔法or欧拉法求解微分⽅程matlab实现举例：分别⽤欧拉法和龙哥库塔法求解下⾯的微分⽅程我们知道的欧拉法（Euler）"思想是⽤先前的差商近似代替倒数"，直⽩⼀些的编程说法即：f(i+1)=f(i)+h*f(x,y)其中h是设定的迭代步长，若精度要求不⾼，⼀般可取0.01。
在定义区间内迭代求解即可。
龙哥库塔法⼀般⽤于⾼精度的求解，即⾼阶精度的改进欧拉法，常⽤的是四阶龙哥库塔，编程语⾔如下：y(i+1)=y(i)+h*(k1+2*K2+2*k3+k4)/6;k1=f(xi,yi)k2=f(xi+h/2,yi+h*k1/2);k3=f(xi+h/2,yi+h*k2/2);k4=f(xi+h,yi+h*k3);设置终⽌条件迭代求解。
matlab实现程序如下：%% 龙哥库塔or欧拉法求解微分⽅程t=0:0.01:3; %⾃变量范围f = [];g = [];f(1) = 0.1; %f初值g(1) = 1; %g初值h=0.001;for i=1:length(t)%% 欧拉法% f(i+1) =f(i)+h*(exp(f(i)*sin(t(i)))+g(i));% g(i+1) =g(i)+h*(exp(g(i)*cos(t(i)))+f(i));%% 龙哥库塔kf1 = exp(f(i)*sin(t(i)))+g(i);%g(i)相当于常值kf2 = exp((f(i)+kf1*h/2)*sin(t(i)+h/2))+g(i);kf3 = exp((f(i)+kf2*h/2)*sin(t(i)+h/2))+g(i);kf4 = exp((f(i)+kf3*h)*sin(t(i)+h))+g(i);f(i+1) = f(i) + h*(kf1+2*kf2+2*kf3+kf4)/6;kg1 = exp(f(i)*cos(t(i)))+f(i);%f(i)相当于常值kg2 = exp((f(i)+kg1*h/2)*cos(t(i)+h/2))+f(i);kg3 = exp((f(i)+kg2*h/2)*cos(t(i)+h/2))+f(i);kg4 = exp((f(i)+kg3*h)*cos(t(i)+h))+f(i);g(i+1) = g(i) + h*(kg1+2*kg2+2*kg3+kg4)/6;endfigure(1)plot(t,f(1:length(t)),t,g(1:length(t)));legend('f函数','g函数')。
Matlab作业龙格库塔欧拉方法解二阶微分方程精选
? 现在求解的是一个类似的问题，在这里的单摆是一种特别 的单摆，具有均匀的质量M分布在长为2的臂状摆上。
? 使用能量法（动能定理）建立方程 T ? W
1 J? 2 ? mg? h
2
? 化简得到
d 2?
dt 2
?
7.35499cos? （重力加速度取9.80665m/s2）
计算
? 边值条件y(0)=0,y'(0)=0.
接下来进行对比
计算
? 运行第三个程序：在一幅图中显示欧拉法和RK4法，随 着截断误差的积累，欧拉法产生了较大的误差
? h=0.01
h=0.0001
总结
? 通过这两种方法计算出角度峰值y=3.141593，周期是 1.777510。
? Euler方法结构简单，但是由于截断误差，使误差较大。 ? RK4是很好的方法，很稳定，由于到五阶的时候精度并
没有相应提升，所以四阶是很常用的方法。
左平成 S14060663 储建研14-2 理论力学专业
Matlab 应用
使用Euler和Rungkutta方法解 臂状摆的能量方程
1. 背景
? 定 理
M ? J?
? 化简得到
d 2?
dt 2
?
g l
sin ?
?
0
? 这化样为在，小这于样5度比的较时容候易容解易。si实n?际上简
这是一?个解二阶常微分方程的问题。
2. 问题
? 1. 使用Euler方法
? 精度随着h的减小而更高，因为向前欧拉方法的整体截 断误差与h同阶，（因为用了泰勒公式）所以欧拉方法 的稳定区域并不大。通过减小h增加了稳定性。
h=0.01
h=0.0001
matlab软件欧拉算法教程
y( xn1 ) y( xn ) hK
*
寻求计算平均斜率的算法
考察欧拉法，以xn的斜率值
K1 f ( xn , yn )
作为平均斜率
考察改进的欧拉法，可以将其改写为：
1 1 yi 1 yn h K1 K 2 2 2 K1 f ( x n , yn ) K 2 f ( xn h, yn hK1 )
xn
h f ( x , y( x ))dx f ( xn , y( xn )) f ( xn1 , y( xn1 )) O( h3 ) 2
h yn1 yn [ f ( xn , yn ) f ( xn1 , yn1 )] 梯形格式 2
梯形格式是显式Euler格式与隐式Euler格式的 算术平均 Euler格式是显式算法，计算量小，但精度低 梯形格式，精度较高，但是隐式算法，需要通过 迭代过程求解，计算量大
四、四阶Runge-Kutta方法
继续上述过程，可以进一步导出四阶Runge-Kutta 格式
yn1 yn h K1 2 K 2 2 K 3 K 4 , 6 K f x , y , n n 1 h K 2 f x n 1 , yn K1 , 2 2 K f x , y h K , n 1 n 2 3 2 2 K 4 f xn1 , yn hK 3 . Biblioteka 利用数值积分求积分项的方法离散
(1)左矩形法
xn1
f ( x , y( x ))dx hf ( xn , y( xn )) O( h2 )
xn
y( xn1 ) y( xn ) hf ( xn , y( xn ))
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)仍然收敛，但是收敛速度大大降低。
解微分方程欧拉法,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决定。
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下的应用
科教论坛 !"#!$%&$'(') *+&,-./&$01$21(3$&)%)$$%%(3
微分方程数值解之欧拉法在 J7<97F下的应用
梁春叶4王桥明4孙远通4叶晓艳4曾宝莹
玉林师范学院数学与统计学院!广西玉林!(#+$$$
摘4要微分方程在实际应用中十分广泛涉及领域众多但对于微分方程的数值解的计算仍然有很大挑战 本文着重 对微分方程数值解求解的常用的一类基础方法欧拉法进行在 EH#[HC的应用下的一个简单介绍
求出$ 以下为在 50;"0)中实现!
F8/F68& 6'
+'
Copyright©博看网 . All Rights Reserved.
科教论坛
%I&Fh(J>>& 8'hh63t3i8326% KG&( h+ 8&$'hh( 27, % 8h(FG*H%& %I&F#KG&(' 由此得该方程的精确解为!
科技风 "#"$ 年 % 月
8hFJ/E*J>8& 8' 由此得解析解为!8h67 i63 N6N$' $
)% $) $) )% 使用欧拉法求解#如下! A h%3$% >hu& 6#8' 63t(i63t)% 6h+$!A!$%, % -hFJ9%& 6#)' i$% 8$ h+) 27#9%$GF&$#-' , % >G$& h$!8$& &N$'h8$& &' NA!>& 6& &' #8$& &' ' % %&( 8h63t732)%i63t332$)N632$)N$' 2)%% E*G7& 6#8#pGBp#6#8$#pN@p' *%L%&(& p,6?K7p#p,+*%$p' 运行程序可输出结果为!
matlab求解微分方程数值解与解析解
matlab求解微分方程数值解与解析解微分方程是数学中的重要内容，它描述了物理、工程、经济等领域中的许多现象和问题。
在实际应用中，我们经常需要求解微分方程的解析解或数值解。
本文将以Matlab为工具，探讨如何求解微分方程并对比解析解与数值解的差异。
一、引言微分方程是描述自然界中许多现象和问题的数学语言，它包含了未知函数及其导数与自变量之间的关系。
微分方程的求解可以帮助我们了解问题的性质和变化规律，并为实际应用提供参考。
在许多情况下，微分方程的解析解很难求得，这时我们可以利用计算机进行数值求解。
二、微分方程的数值解法1.欧拉法欧拉法是最简单的数值求解微分方程的方法之一。
它通过将微分方程转化为差分方程，然后利用离散的点逼近连续的解。
具体步骤如下：（1）将微分方程转化为差分方程，即用近似的导数代替真实的导数；（2）选择初始条件，即确定初始点的值；（3）选择步长和求解区间，即确定求解的范围和步长；（4）使用迭代公式计算下一个点的值；（5）重复步骤（4），直到达到指定的求解区间。
2.改进的欧拉法欧拉法存在精度较低的问题，为了提高精度，可以使用改进的欧拉法。
改进的欧拉法是通过使用两次导数的平均值来计算下一个点的值，从而提高了数值解的精度。
3.龙格-库塔法龙格-库塔法是一种常用的数值求解微分方程的方法，它通过使用多个点的导数来逼近连续解。
龙格-库塔法的步骤如下：（1）选择初始条件和步长；（2）使用迭代公式计算下一个点的值；（3）计算下一个点的导数；（4）根据导数的值和步长计算下一个点的值；（5）重复步骤（3）和（4），直到达到指定的求解区间。
龙格-库塔法的精度较高，适用于求解一阶和高阶微分方程。
三、微分方程的解析解解析解是指能够用公式或函数表示的方程的解。
有些微分方程具有解析解，可以通过数学方法求得。
例如，一阶线性常微分方程和某些特殊类型的二阶微分方程等。
解析解的优势在于精确性和直观性，能够帮助我们深入理解问题的本质。
7-17-2微分方程01欧拉法02R-K法
1 y ′′ (ξ ) h 2 2
二阶泰勒公式
y ( xn +1 ) = y ( xn ) + y ′ ( xn ) h +
两式相减， 两式相减，由设 yn=y(xn ) ，有
h2 y ( xn +1 ) − yn +1 = y′′ (ξ ) = O ( h 2 ) 2
一阶常微分方程初 值问题
dy = f ( x, y ) dx y ( x0 ) = y0
的数值解法。 基 本思想 ：常 微分方 程初值 问题的 数值解 是求 微分方程 的 解 y ( x) 即 微 分 方 程 初 值 问 题 的 积 分 曲 线 ） 在 区 间 [ a , b] 中 （ ， 点（ 给 定 一 系 列 点 （ 节 点 ） x n = x n − 1 + hn （ n = 1, 2, L ） 上 的 近 似 值 y n 。 这 里 hn 为 x n −1 到 x n 的 步 长 ， 且 h n > 0 。
等距， 设 x0 , x1 , x2 ,⋅ ⋅ ⋅, xn 等距，步长
y′( x ) ≈
y ( x + h ) ≈ y ( x ) + h y ′( x ) = y ( x ) + h f ( x , y )
令x=xn , x+h=xn+1 , y(xn)≈yn ，y(xn+1 ) ≈yn+1 ,
初值问题离散化为
h2 y ( xn +1 ) − yn +1 ≈ − y ''( xn ) 2
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一、什么是常微分方程初值问题常微分方程是研究函数、导数和变量之间的关系的方程，初值问题则是给定函数在某一点的初始值，通过求解常微分方程可以得到函数在其他点的数值解。
常微分方程初值问题作为数值计算领域的重要问题，对于很多实际应用具有重要意义。
二、欧拉法简介欧拉法是求解常微分方程初值问题的一种数值计算方法之一。
它的基本思想是通过离散化函数的导数，将微分方程转化为差分方程，然后通过迭代计算逼近函数在其他点的近似解。
欧拉法虽然简单，但在一定条件下可以得到较好的近似解。
三、使用M atlab实现常微分方程初值问题的欧拉法在M at la b中，我们可以使用以下步骤来实现常微分方程初值问题的欧拉法：1.定义微分方程首先，需要定义待求解的常微分方程。
我们可以使用匿名函数来定义方程，例如：f=@(t,y)y-t^2+1;这里，`t`表示自变量，`y`表示因变量。
2.定义初值和步长接下来，需要定义初始条件和步长。
初始条件指定了函数在某一点的初始值，步长表示要计算的点之间的间隔。
例如：t0=0;%初始时间y0=0.5;%初始值h=0.1;%步长3.进行迭代计算使用欧拉法进行迭代计算的核心步骤如下：f o ri=1:Nt(i+1)=t(i)+h;%计算下一个时间点y(i+1)=y(i)+h*f(t(i),y(i));%计算下一个函数值e n d其中，`N`是迭代次数，`t`和`y`是存储计算结果的数组。
4.绘制图像最后，我们可以使用M at la b的绘图函数将结果可视化，例如：p l ot(t,y,'o-');x l ab el('t');y l ab el('y');t i tl e('常微分方程初值问题的欧拉法解');四、总结通过使用M at la b，我们可以方便地实现常微分方程初值问题的欧拉法。
首先需要定义微分方程、初值和步长，然后通过迭代计算可以得到函数在其他点的数值解。
matlab用欧拉法求常微分方程初值
Matlab中欧拉法求解常微分方程初值问题一、概念介绍在数学和工程领域，常微分方程初值问题是一个广泛应用的数学概念。
它描述了一个未知函数在给定初始条件下的行为。
而欧拉法则是一种常用的数值方法，用来解决常微分方程初值问题。
在Matlab中，我们可以利用欧拉法来求解常微分方程问题，从而得到函数在给定初始条件下的近似解。
二、欧拉法的基本原理欧拉法的基本思想是通过离散化微分方程，将其转化为递推的差分方程。
考虑一个一阶常微分方程初值问题：\[ \frac{dy}{dx} = f(x, y), \quad y(x_0) = y_0 \]在欧拉法中，我们采用递推的方式，根据已知的初始条件和微分方程的性质，通过迭代来得到逼近解的数值结果。
具体地，我们首先将自变量$x$的范围进行等间距分割，得到$x_0, x_1, x_2, ..., x_n$，并将步长记为$h$。
根据微分方程的性质，我们可以根据已知的初始条件$y(x_0) = y_0$，通过迭代计算得到近似解$y(x_1), y(x_2), ..., y(x_n)$。
三、Matlab中的欧拉法求解在Matlab中，我们可以利用欧拉法来求解常微分方程初值问题。
以求解一阶常微分方程为例，假设我们需要求解以下的常微分方程初值问题：\[ \frac{dy}{dx} = -2xy, \quad y(0) = 1 \]我们可以利用欧拉法的思想，将自变量$x$的范围进行离散化，然后根据欧拉法的递推公式，利用迭代的方式得到近似解的数值结果。
具体地，在Matlab中，我们可以编写如下代码来实现欧拉法的求解过程：```matlabfunction y = euler_method(f, x0, y0, h, n)% 初始化存储结果的数组x = zeros(1, n+1);y = zeros(1, n+1);% 将初始条件存入数组x(1) = x0;y(1) = y0;% 利用欧拉法进行迭代for i = 1:nx(i+1) = x(i) + h;y(i+1) = y(i) + h * f(x(i), y(i));end% 返回近似解的数值结果plot(x, y); % 绘制解的图像end```在上述代码中，我们定义了一个名为`euler_method`的函数，其中包含了欧拉法的计算过程。
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 Euler法解常微分方程
Euler 法解常微分方程Euler 法解常微分方程算法：Step 1 分别取积分上限、积分下限、步长Step 2计算h n n +=判断b n ≤是否成立，成立转到Step 3，否则继续进行Step 4 Step 3 计算),(1n n n n y x hf y y +=+Step 4 ),(1n n n n y x hf y y +=+Euler 法解常微分方程算程序：function euler2(fun,y0,A,h)%fun--y'%y0---初值%A----x 取值范围%a----x 左区间端点值%b----x 右区间端点值%h----给定步长x=min(A);b=max(A);y=y0;while x<b-hb=y;y=y+h*feval(fun,x,b)x=x+h;end例：用Euler 法计算下列初值问题（取步长h=0.2）)6.00(1)0('2≤≤⎩⎨⎧=--=x y xy y y 输入：fun=inline('-y-x*y^2')euler2(fun,1,[0 0.6],0.2)得到：y =0.8000y =0.6144y =0.4613指导教师： 年 月 日改进Euelr 法解常微分方程改进Euler 法解常微分方程算法：Step 1 分别取积分上限、积分下限、步长Step 2 取一个以h 为步长，a ，b 分别为左右端点的矩阵Step 3 （1）做显性Euler 预测),(1n n i i y x hf y y +=+（2）将1+i y 带入(),([2h 11++++=i i i i i x f y x f y y Step 4计算h n n +=判断b n ≤是否成立，成立返回Step 5 )],(),([2h 111+++++=i i i i i i y x f y x f y y 改进Euler 法解常微分方程算程序：function gaijineuler2(fun,y0,A,h)%fun--y'%y0---初值%A----x 取值范围%a----x 左区间端点值%b----x 右区间端点值%h----给定步长a=min(A);b=max(A);x=a:h:b;y(1)=y0;for i=1:length(x)-1w1=feval(fun,x(i),y(i));y(i+1)=y(i)+h*w1;w2=feval(fun,x(i+1),y(i+1));y(i+1)=y(i)+h*(w1+w2)/2;endx=x'y=y'例：用改进Euler 法计算下列初值问题（取步长h=0.25）)50(2)0('2≤≤⎩⎨⎧=-=x y xy y 输入：fun=inline('-x*y^2')gaijineuler2(fun,2,[0 5],0.25)得到：x =0.25000.50000.75001.00001.25001.50001.75002.00002.25002.50002.75003.00003.25003.50003.75004.00004.25004.50004.75005.0000y =2.00001.87501.59391.28241.00960.79320.62820.50370.40970.33790.28240.23890.20430.17650.15380.13520.11960.10660.09550.08610.0779指导教师：年月日。
matlab-常微分方程数值解法
第二十页，编辑于星期一：二十三点 三十分。
2.2 改进的欧拉近似方法
对于一般的常微分方程：
dy
d
x
f (x, y)
y ( x 0 ) y 0
节点： x0,x1, ,xi, ,xn, 步长： hxi1xi
例如,质点的加速运动，简谐振动等。
F m dv dt
d2x 2x2 0
dt 2
简单问题可以求得解析解，多数实际问题靠数值求解。
第四页，编辑于星期一：二十三点 三十分。
一阶常微分方程(ODE )初值问题 ：
ODE ：Ordinary Differential Equation
dy
f
(x,
y)
dx
y0
y1 y2 …. yn ….
第八页，编辑于星期一：二十三点 三十分。
在xn节点上，微分方程可以写为
y ( x n 1 ) y ( x n ) f y ( x n ) ,x n h
作如下近似：
yn y(tn)
则得到欧拉解法递推公式的一般形式：
y n 1 y n f(y n,x n )h
② 用向后的欧拉方法校正---迭代
y n (k 1 1 ) y n h f(x n 1 ,y n (k 1 ))
第二十三页，编辑于星期一：二十三点 三十分。
2.2.2 改进的欧拉方法
积分法求欧拉公式时，积分采用梯形近似，即得到改进的欧 拉方法。
y(xn1)y(xn)xxnn1 f(x,y)dx
2
)
k 4 f ( x n h , y n h k 3 )
数值分析Matlab作业龙格库塔欧拉方法解二阶微分方程
Matlab 应用使用Euler 和Rungkutta 方法解臂状摆的能量方程背景 单摆是常见的物理模型，为了得到摆角θ的关于时间的函数，来描述单摆运动。
由角动量定理我们知道εJ M =化简得到0sin 22=+θθl g dt d 在一般的应用和计算中，只考虑摆角在5度以内的小摆动，因为可以吧简化为，这样比较容易解。
实际上这是一个解二阶常微分方程的问题。
θ在这里的单摆是一种特别的单摆，具有均匀的质量M 分布在长为2的臂状摆上，使用能量法建立方程W T =h mg ∆=2J 21ω化简得到四阶龙格库塔方法使用四级四阶经典显式Rungkutta公式稳定性很好，RK4法是四阶方法，每步的误差是h5阶，而总积累误差为h4阶。
所以比欧拉稳定。
运行第三个程序：在一幅图中显示欧拉法和RK4法，随着截断误差的积累，欧拉法产生了较大的误差h=0.01h=0.0001，仍然是开始较为稳定，逐渐误差变大总结：RK4是很好的方法，很稳定，而且四阶是很常用的方法，因为到五阶的时候精度并没有相应提升。
通过这两种方法计算出角度峰值y=3.141593，周期是1.777510。
三个程序欧拉法clear;clch=0.00001;a=0;b=25;x=a:h:b;y(1)=0;z(1)=0;for i=1:length(x)-1 % 欧拉y(i+1)=y(i)+h*z(i);z(i+1)=z(i)+h*7.35499*cos(y(i));endplot(x,y,'r*');xlabel('时间');ylabel('角度');A=[x,y];%y(find(y==max(y)))%Num=(find(y==max(y)))[y,T]=max(y);fprintf('角度峰值等于%d',y) %角度的峰值也就是πfprintf('\n')fprintf('周期等于%d',T*h)%周期legend('欧拉');龙格库塔方法先定义函数rightf_sys1.mfunction w=rightf_sys1(x,y,z)w=7.35499*cos(y);clear;clc;%set(0,'RecursionLimit',500)h=0.01;a=0;b=25;x=a:h:b;RK_y(1)=0; %初值%RK_z(1)=0;初值for i=1:length(x)-1K1=RK_z(i); L1=rightf_sys1(x(i),RK_y(i),RK_z(i));%K1 and L1K2=RK_z(i)+0.5*h*L1;L2=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K1,RK_z(i)+0.5*h*L1);K3=RK_z(i)+0.5*h*L2;L3=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K2,RK_z(i)+0.5*h*L2);K4=RK_z(i)+h*L3;% K4L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3);and L4RK_y(i+1)=RK_y(i)+1/6*h*(K1+2*K2+2*K3+K4);RK_z(i+1)=RK_z(i)+1/6*h*(L1+2*L2+2*L3+L4);endplot(x,RK_y,'b+');xlabel('Variable x');ylabel('Variable y');A=[x,RK_y];[y,T]=max(RK_y);legend('RK4方法');fprintf('角度峰值等于%d',y) %角度的峰值也就是πfprintf('\n')%周期fprintf('周期等于%d',T*h)两个方法在一起对比使用跟上一个相同的函数rightf_sys1.mclear;clc; %清屏h=0.0001;a=0;b=25;x=a:h:b;Euler_y(1)=0;%欧拉的初值Euler_z(1)=0;RK_y(1)=0;%龙格库塔初值RK_z(1)=0;for i=1:length(x)-1%先是欧拉法Euler_y(i+1)=Euler_y(i)+h*Euler_z(i);Euler_z(i+1)=Euler_z(i)+h*7.35499*cos(Euler_y(i));%龙格库塔K1=RK_z(i); L1=rightf_sys1(x(i),RK_y(i),RK_z(i)); % K1 andL1K2=RK_z(i)+0.5*h*L1;L2=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K1,RK_z(i)+0.5*h*L1);% K2 and L2K3=RK_z(i)+0.5*h*L2;L3=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K2,RK_z(i)+0.5*h*L2);% K3 and L3K4=RK_z(i)+h*L3; L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3); K4 and L4RK_y(i+1)=RK_y(i)+1/6*h*(K1+2*K2+2*K3+K4);RK_z(i+1)=RK_z(i)+1/6*h*(L1+2*L2+2*L3+L4);endplot(x,Euler_y,'r-',x,RK_y,'b-');[y,T]=max(RK_y);%角度的峰值也就是πfprintf('角度峰值等于%d',y)fprintf('\n')%周期fprintf('周期等于%d',T*h)xlabel('时间');ylabel('角度');legend('欧拉','RK4');。
Matlab实验报告五(微分方程求解Euler折线法)
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))
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
欧拉方法（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决定。
当四个斜率取平均时，中点的斜率有更大的权值：
RK4法是四阶方法，也就是说每步的误差是h5阶，而总积累误差为h4阶。
注意上述公式对于标量或者向量函数(y可以是向量)都适用。
例子：
下面给出了数值求解该微分方程的简单程序。
其中y1,y2,y3,y4分别为向前欧拉公式，改进的欧拉公式，4级4阶龙格-库塔公式及精确解。
h=;
x=0:h:1;
y1=zeros(size(x));
y1(1)=1;
y2=zeros(size(x));
y2(1)=1;
y3=zeros(size(x));
y3(1)=1;
for i1=2:length(x)
y1(i1)=y1(i1-1)+h*(y1(i1-1)-2*x(i1-1)/y1(i1-1));
k1=y2(i1-1)-2*x(i1-1)/y2(i1-1);
k2=y2(i1-1)+h*k1-2*x(i1)/(y2(i1-1)+h*k1);
y2(i1)=y2(i1-1)+h*(k1+k2)/2;
k1=y2(i1-1)-2*x(i1-1)/y2(i1-1);
k2=y2(i1-1)+h*k1/2-2*(x(i1-1)+h/2)/(y2(i1-1)+h*k1/2);
k3=y2(i1-1)+h*k2/2-2*(x(i1-1)+h/2)/(y2(i1-1)+h*k2/2);
k4=y2(i1-1)+h*k3-2*(x(i1-1)+h)/(y2(i1-1)+h*k3);
y3(i1)=y3(i1-1)+(k1+2*k2+2*k3+k4)*h/6;
end
y4=sqrt(1+2*x);
%plot(x,y1,x,y2,x,y3,x,y4)
%legend('y1','y2','y3','y4')
plot(x,y4-y1,x,y4-y2,x,y4-y3)
legend('y1','y2','y3')。