MATLAB改进欧拉法与四阶龙格-库塔求解一阶常微分方程
MATLAB改进欧拉法与四阶龙格-库塔求解一阶常微分方程
姓名：樊元君学号：02 日期：一、实验目的掌握MATLAB语言、C/C++语言编写计算程序的方法、掌握改进欧拉法与四阶龙格-库塔求解一阶常微分方程的初值问题。
掌握使用MATLAB程序求解常微分方程问题的方法。
：二、实验内容1、分别写出改进欧拉法与四阶龙格-库塔求解的算法，编写程序上机调试出结果，要求所编程序适用于任何一阶常微分方程的数值解问题，即能解决这一类问题，而不是某一个问题。
实验中以下列数据验证程序的正确性。
求，步长h=。
*2、实验注意事项的精确解为，通过调整步长，观察结果的精度的变化^)三、程序流程图：●改进欧拉格式流程图：~|●四阶龙格库塔流程图：]四、源程序：●改进后欧拉格式程序源代码：function [] = GJOL(h,x0,y0,X,Y)format longh=input('h=');…x0=input('x0=');y0=input('y0=');disp('输入的范围是：');X=input('X=');Y=input('Y=');n=round((Y-X)/h);\i=1;x1=0;yp=0;yc=0;for i=1:1:nx1=x0+h;yp=y0+h*(-x0*(y0)^2);%yp=y0+h*(y0-2*x0/y0);%·yc=y0+h*(-x1*(yp)^2);%yc=y0+h*(yp-2*x1/yp);%y1=(yp+yc)/2;x0=x1;y0=y1;y=2/(1+x0^2);%y=sqrt(1+2*x0);%fprintf('结果=%.3f,%.8f,%.8f\n',x1,y1,y);:endend●四阶龙格库塔程序源代码：function [] = LGKT(h,x0,y0,X,Y)。
format longh=input('h=');x0=input('x0=');y0=input('y0=');disp('输入的范围是：');"X=input('X=');Y=input('Y=');n=round((Y-X)/h);i=1;x1=0;k1=0;k2=0;k3=0;k4=0;for i=1:1:n~x1=x0+h;k1=-x0*y0^2;%k1=y0-2*x0/y0;%k2=(-(x0+h/2)*(y0+h/2*k1)^2);%k2=(y0+h/2*k1)-2*(x0+h/2)/(y0+h/2*k1);% k3=(-(x0+h/2)*(y0+h/2*k2)^2);%k3=(y0+h/2*k2)-2*(x0+h/2)/(y0+h/2*k2);% k4=(-(x1)*(y0+h*k3)^2);%k4=(y0+h*k3)-2*(x1)/(y0+h*k3);%…y1=y0+h/6*(k1+2*k2+2*k3+k4);%y1=y0+h/6*(k1+2*k2+2*k3+k4);%x0=x1;y0=y1;y=2/(1+x0^2);%y=sqrt(1+2*x0);%fprintf('结果=%.3f,%.7f,%.7f\n',x1,y1,y);end·end*五、运行结果：改进欧拉格式结果：;}四阶龙格库塔结果：步长分别为：和时，不同结果显示验证了步长减少，对于精度的提高起到很大作用，有效数字位数明显增加。
数值分析常微分方程求解实验
实验报告
一、实验目的
解初值问题各种方法比较。
二、实验题目
给定初值问题
⎪⎩⎪⎨⎧=≤<+=,
0)1(,21y x xe x y dx dy x ， 其精确解为)(e e x y x -=，按
（1）改进欧拉法，步长01.0,1.0==h h ；
（2）四阶标准龙格-库塔法，步长1.0=h ；
求在节点)10,...,2,1(1.01=+=k k x k 处的数值解及其误差，比较各个方法的优缺点。
三、实验原理
改进欧拉法程序，四阶标准龙格-库塔法程序。
四、实验内容及结果
五、实验结果分析
实验2中改进欧拉法和四阶标准龙格-库塔法的比较：
结果的第一个ans是x的值与对应的y的值，第二个ans是精确解x的对应值y，第三个ans 是与精确值的误差百分数。
通过误差百分数的比较，可以明显的发现改进欧拉法比四阶标准龙格-库塔法更精确。
微分方程的数值解法matlab(四阶龙格—库塔法)
解析解： x x x1 3 2(((ttt))) 0 .0 8 1 1 2 P k 8 0siw n t) (2 .6 3 0 3 3 P k 0siw n t) (0 .2 12 2 2 P k 0siw n t)(
第一个质量的位移响应时程
Y (t)A(Y t)P(t)
(2)
Y (t)A(Y t)P(t)
3. Matlab 程序(主程序：ZCX)
t0;Y0;h;N;P0,w; %输入初始值、步长、迭代次数、初始激励力；
for i = 1 : N
t1 = t0 + h
P=[P0*sin(w*t0);0.0;0.0]
%输入t0时刻的外部激励力
Van der Pol方程
% 子程序 （程序名： dYdt.m ） function Ydot = dYdt (t, Y) Ydot=[Y(2);-Y(2)*(Y(1)^2-1)-Y(1)];
或写为
function Ydot = dYdt (t, Y) Ydot=zeros(size(Y)); Ydot(1)=Y(2); Ydot(2)=-Y(2)*(Y(1).^2-1)-Y(1)];
Solver解算指令的使用格式
说明：
t0：初始时刻；tN：终点时刻 Y0：初值； tol：计算精度
[t, Y]=solver (‘ODE函数文件名’, t0, tN, Y0, tol);
ode45
输出宗量形式
y1 (t0 )
Y
y1
(t1
)
y
1
(t
2
)
y2 (t0 )
y
2
(
t1
)
y
2
(
t
欧拉法求解一阶微分方程matlab
为了更好地理解欧拉法求解一阶微分方程在Matlab中的应用，我们首先来了解一些背景知识。
一阶微分方程是指只含有一阶导数的方程，通常表示为dy/dx=f(x,y)，其中f(x,y)是关于x和y的函数。
欧拉法是一种常见的数值解法，用于求解微分方程的近似数值解。
它是一种基本的显式数值积分方法，通过将微分方程转化为差分方程来进行逼近。
在Matlab中，我们可以利用欧拉法求解一阶微分方程。
我们需要定义微分方程的函数表达式，然后选择合适的步长和初始条件，最后使用循环计算逼近解。
下面我们来具体讨论如何在Matlab中使用欧拉法来求解一阶微分方程。
我们假设要求解的微分方程为dy/dx=-2x+y，初始条件为y(0)=1。
我们可以通过以下步骤来实现：1. 我们需要在Matlab中定义微分方程的函数表达式。
在Matlab中，我们可以使用function关键字来定义函数。
在这个例子中，我们可以定义一个名为diff_eqn的函数，表示微分方程的右侧表达式。
在Matlab中，这个函数可以定义为：```matlabfunction dydx = diff_eqn(x, y)dydx = -2*x + y;end```2. 我们需要选择合适的步长和初始条件。
在欧拉法中，步长的选择对于数值解的精度非常重要。
通常情况下，可以先尝试较小的步长，然后根据需要进行调整。
在这个例子中，我们可以选择步长h=0.1，并设置初始条件x0=0，y0=1。
3. 接下来，我们可以使用循环来逼近微分方程的数值解。
在每一步，根据欧拉法的迭代公式y(i+1) = y(i) + h * f(x(i), y(i))，我们可以按照下面的Matlab代码计算逼近解：```matlabh = 0.1; % 步长x = 0:h:2; % 定义计算区间y = zeros(1, length(x)); % 初始化y的值y(1) = 1; % 设置初始条件for i = 1:(length(x)-1) % 欧拉法迭代y(i+1) = y(i) + h * diff_eqn(x(i), y(i));end```通过上述步骤，在Matlab中就可以用欧拉法求解一阶微分方程。
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。
4阶RungeKutta法求解一阶常微分方程
《MATLAB语言及应用》大作业姓名:学号:学院:班级:题目编号:2013年10月134阶Runge-Kutta法求解一阶常微分方程。
一、Runge-Kutta法的数学理论龙格-库塔(Runge-Kutta)方法就是一种在工程上应用广泛的高精度单步算法。
由于此算法精度高,采取措施对误差进行抑制,所以其实现原理也较复杂。
该算法就是构建在数学支持的基础之上的。
龙格库塔方法的理论基础来源于泰勒公式与使用斜率近似表达微分,它在积分区间多预计算出几个点的斜率,然后进行加权平均,用做下一点的依据,从而构造出了精度更高的数值积分计算方法。
如果预先求两个点的斜率就就是二阶龙格库塔法,如果预先取四个点就就是四阶龙格库塔法。
一阶常微分方程可以写作:y'=f(x,y),使用差分概念。
(Yn+1-Yn)/h= f(Xn,Yn)推出(近似等于,极限为Yn')Yn+1=Yn+h*f(Xn,Yn)另外根据微分中值定理,存在0<t<1,使得Yn+1=Yn+h*f(Xn+th,Y(Xn+th))这里K＝f(Xn+th,Y(Xn+th))称为平均斜率,龙格库塔方法就就是求得K的一种算法。
利用这样的原理,经过复杂的数学推导(过于繁琐省略),可以得出截断误差为O(h^5)的四阶龙格库塔公式:K1＝f(Xn,Yn);K2=f(Xn+h/2,Yn+(h/2)*K1);K3=f(Xn+h/2,Yn+(h/2)*K2);K4=f(Xn+h,Yn+h*K3);Yn+1=Yn+h*(K1+2K2+2K3+K4)*(1/6)二、Runge-Kutta的算法与流程图在龙格-库塔法中,四阶龙格-库塔法的局部截断误差约为0(h5),被广泛应用于解微分方程的初值问题。
其算法公式为:[]()()()()112341213243/622,/2,/2*,/2,/2*,,n n n n n n n n n n y y h k k k k k f x y k f x h y h k k f x h y h k k f x h y hk +⎧⎫=++++⎪⎪=⎪⎪⎪⎪=++⎨⎬⎪⎪=++⎪⎪⎪⎪=++⎩⎭流程图:(1)、四阶龙格-库塔方法流程图:(2)、实例求解流程图:三、Runge-Kutta的Matlab实现function [x,y]=runge_kutta1(ufunc,y0,h,a,b)%参数表顺序依次就是微分方程组的函数名称,初始值向量,步长,时间起点,时间终点(参数形式参考了ode45函数)n=floor((b-a)/h);%求步数x(1)=a;%时间起点y(:,1)=y0;%赋初值,可以就是向量,但就是要注意维数for ii=1:nx(ii+1)=x(ii)+h;k1=ufunc(x(ii),y(:,ii));k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2);k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2);k4=ufunc(x(ii)+h,y(:,ii)+h*k3);y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6;%按照龙格库塔方法进行数值求解end四、Runge-Kutta的算例实现例:求解常微分方程d(y)/d(x)=-2*y+2*x*x+2*x ,0≤x≤0、5 ,y(0)=1、编辑m文件Untitled3:fun=inline('-2*y+2*x*x+2*x');[x,y]=ode45(fun,[0,0、5],1)>> Untitled3x =0、01250、02500、03750、0625 0、0750 0、0875 0、1000 0、1125 0、1250 0、1375 0、1500 0、1625 0、1750 0、1875 0、2000 0、2125 0、2250 0、2375 0、2500 0、2625 0、2750 0、2875 0、3000 0、31250、33750、35000、36250、37500、38750、40000、41250、42500、43750、45000、46250、47500、48750、5000 y =1、00000、97550、95190、9073 0、8864 0、8663 0、8471 0、8287 0、8112 0、7944 0、7785 0、7633 0、7489 0、7353 0、7224 0、7103 0、6989 0、6883 0、6783 0、6690 0、6605 0、6526 0、6454 0、63880、6277 0、6231 0、6191 0、6157 0、6130 0、6109 0、6093 0、6084 0、6080 0、6083 0、6091 0、6104 0、6124 0、6148 0、6179。
常微分方程的数值解法(欧拉法、改进欧拉法、泰勒方法和龙格库塔法)
［例1］用欧拉方法与改进的欧拉方法求初值问题h 的数值解。
在区间[0,1]上取0.1［解］欧拉方法的计算公式为x0=0;y0=1;x(1)=0.1;y(1)=y0+0.1*2*x0/(3*y0^2);for n=1:9x(n+1)=0.1*(n+1);y(n+1)=y(n)+0.1*2*x(n)/(3*y(n)^2);end;xy结果为x =Columns 1 through 80.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 Columns 9 through 100.9000 1.0000y =Columns 1 through 81.0000 1.0067 1.0198 1.0391 1.0638 1.0932 1.1267 1.1634 Columns 9 through 101.2028 1.2443改进的欧拉方法其计算公式为本题的精确解为()y x=x0=0;y0=1;ya(1)=y0+0.1*2*x0/(3*y0^2);y(1)=y0+0.05*(2*x0/(3*y0^2)+2*x0/(3*ya^2));for n=1:9x(n+1)=0.1*(n+1);ya(n+1)=ya(n)+0.1*2*x(n)/(3*ya(n)^2);y(n+1)=y(n)+0.05*(2*x(n)/(3*y(n)^2)+2*x(n+1)/(3*ya(n+1)^2));end;xy结果为x =Columns 1 through 80.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 Columns 9 through 100.9000 1.0000y =Columns 1 through 81.0000 1.0099 1.0261 1.0479 1.0748 1.1059 1.1407 1.1783 Columns 9 through 101.2183 1.2600［例2］用泰勒方法解x＝0.1, 0.2, …, 1.0处的数值解，并与精确解进行比较。
matlab编写龙格库塔法或欧拉法求解常微分方程数值解
龙格库塔法（Runge-Kutta method）和欧拉法（Euler's method）是两种常用的数值求解常微分方程的方法。
这里分别给出它们的MATLAB实现：1. 龙格库塔法（Runge-Kutta method）：```matlabfunction [y, t] = runge_kutta(f, y0, t0, tf, h)% f: 微分方程函数，输入为[y, t]，输出为dy/dt% y0: 初始值% t0: 初始时间% tf: 结束时间% h: 步长N = round((tf - t0) / h); % 计算迭代次数t = zeros(1, N + 1); % 初始化时间向量y = zeros(size(y0), N + 1); % 初始化解向量t(1) = t0;y(:, 1) = y0;for i = 1:Nk1 = h * f(y(:, i), t(i));k2 = h * f(y(:, i) + k1 / 2, t(i) + h / 2);k3 = h * f(y(:, i) + k2 / 2, t(i) + h / 2);k4 = h * f(y(:, i) + k3, t(i) + h);y(:, i + 1) = y(:, i) + (k1 + 2 * k2 + 2 * k3 + k4) / 6;t(i + 1) = t(i) + h;endend```2. 欧拉法（Euler's method）：```matlabfunction [y, t] = euler_method(f, y0, t0, tf, h)% f: 微分方程函数，输入为[y, t]，输出为dy/dt% y0: 初始值% t0: 初始时间% tf: 结束时间% h: 步长N = round((tf - t0) / h); % 计算迭代次数t = zeros(1, N + 1); % 初始化时间向量y = zeros(size(y0), N + 1); % 初始化解向量t(1) = t0;y(:, 1) = y0;for i = 1:Ny(:, i + 1) = y(:, i) + h * f(y(:, i), t(i));t(i + 1) = t(i) + h;endend```使用这两个函数时，需要定义一个表示微分方程的函数`f`，例如：```matlabfunction dydt = my_ode(y, t)dydt = -y; % 一个简单的一阶线性微分方程：dy/dt = -yend```然后调用相应的求解函数，例如：```matlaby0 = 1; % 初始值t0 = 0; % 初始时间tf = 5; % 结束时间h = 0.1; % 步长[y_rk, t_rk] = runge_kutta(@my_ode, y0, t0, tf, h);[y_euler, t_euler] = euler_method(@my_ode, y0, t0, tf, h);```。
matlab四阶龙格库塔法解方程组
matlab四阶龙格库塔法解方程组【原创实用版】目录1.MATLAB 简介2.四阶龙格 - 库塔法简介3.用 MATLAB 实现四阶龙格 - 库塔法解方程组的步骤4.结论正文1.MATLAB 简介MATLAB 是一种广泛使用的数学软件，它提供了强大的数值计算和数据分析功能。
MATLAB 中有许多现成的函数和工具箱，可以方便地解决各种数学问题。
在工程、科学和金融领域等领域，MATLAB 都有着广泛的应用。
2.四阶龙格 - 库塔法简介四阶龙格 - 库塔法（RK4）是一种常用的数值积分方法，可以用于求解常微分方程组。
该方法具有较高的精度和稳定性，通常比其他低阶方法需要更少的计算步骤。
四阶龙格 - 库塔法的基本思想是将求解过程分为几个步骤，通过对各阶导数进行适当的组合和积分，最终得到方程组的解。
3.用 MATLAB 实现四阶龙格 - 库塔法解方程组的步骤下面是一个简单的示例，展示如何使用 MATLAB 实现四阶龙格 - 库塔法解方程组。
假设我们要求解如下常微分方程组：y" = x^2 + yz" = x + y我们可以按照以下步骤进行：(1) 创建一个 MATLAB 脚本，定义方程组和初始条件。
例如：```matlabfunction dXdt = rk4(t, X, params)% 设置参数h = 0.01; % 时间步长n = 100; % 时间步数X0 = [1; 0; 0]; % 初始条件% 计算 k1, k2, k3, k4k1 = h*(params(1) + params(3));k2 = h*(params(2) + params(4));k3 = h*(params(2) + params(4));k4 = h*(params(1) + params(3));% 循环求解for i = 1:ndXdt = [k1*X(i, 1); k1*X(i, 2); k1*X(i, 3)];X(i+1, :) = X(i, :) + dXdt;dXdt = [k2*X(i+1, 1); k2*X(i+1, 2); k2*X(i+1, 3)]; X(i+1, :) = X(i+1, :) + dXdt;dXdt = [k3*X(i+1, 1); k3*X(i+1, 2); k3*X(i+1, 3)];X(i+1, :) = X(i+1, :) + dXdt;dXdt = [k4*X(i+1, 1); k4*X(i+1, 2); k4*X(i+1, 3)];X(i+1, :) = X(i+1, :) + dXdt;endend```(2) 定义求解函数，并设置时间范围、时间步长等参数。
龙格-库塔法,求解常微分方程
隆格库塔法求解常微分方程摘要科学技术中常常需要求解常微分方程的定解问题，这里问题最简单的形式，是本章将着重考察的一阶方程的初值问题.虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问题中归结出来的微分方程主要靠数值解法求解.本文着重讨论了隆格库塔法求解一阶常微分方程的初值问题，采用了精度较高的经典的四阶隆格库塔法,然后通过对实例运用Matlab编程进行计算求解，为了体现计算结果的精确性和方法的优越性，再采用了欧拉法和预估较正法对实例进行计算求解作为比较.通过比较三种方法的计算精度，发现四阶经典龙格-库塔方法的误差最小，预估较正法其次，欧拉方法误差则比较大.最后通过选取不同的步长，研究了不同的步长对隆格库塔法求解常微分方程初值问题的计算精度的影响.总之，本文全面分析了隆格库塔法在求解常微分方程的应用，相比与其他的数值解法，隆格库塔法计算精度较高，收敛性较好，其中四阶的隆格库塔法的效率最高，精度也最高.关键词：四阶隆格库塔法；欧拉法；预估较正法；一阶常微分方程；MatlabRunge Kutta Method For Solving Ordinary Differential EquationsABSTRACTProblem solving ordinary differential equations are often needed in science andtechnology. the problem in the simplest form is the initial value problem of first order equations in this chapter ,which will be discussed. Although there are various analytical methods for solving ordinary differential equations, the analytical method can only be used to solve some special types of equations.differential equations can be summed up the actual problems whichThis paper discusses the initial value problem of Runge Kutta Barclays by solving a differential equation, using the four order Runge Kutta method with high accuracy.for instance through classic Matlab programming calculation, the superiority in order to accurately and reflect the calculation result, then the Euler method and the prediction correction method for instance by calculation through the calculation precision. The comparison of three kinds of methods, found that the error of four order Runge Kutta method of minimum, prediction correction method secondly, Euler method error is relatively large. Finally, by selecting different step, study the affect the calculation accuracy of different step of Runge Kutta method to solve initial value problems of ordinary differential equations.In short, this paper comprehensively analyzes the application of Runge Kutta method for solving ordinary differential equations, compared with the numerical solution of other, higher accuracy Runge Kutta method, good convergence, the Runge Kutta method of order four of the highest efficiency and its precision is the highest.Key words: Four order Runge Kutta method; Euler method; prediction correction method; first order ordinary differential equations; Matlab目录1 问题的提出 (1)1.1 问题背景............................................... . (1)1.2 问题的具体内容 (1)2 问题假设 (2)3 符号系统 (2)4 问题的分析 (3)4.1 欧拉格式 (3)4.2 预估较正法 (3)4.3 四阶隆格库塔法的格式 (4)5 模型的建立与求解 (4)5.1 隆格库塔法的基本原理 (4)5.1.1 Taylor级数 (4)5.1.2 隆格库塔法的基本思想 (4)5.1.3 四阶的隆格库塔法 (5)5.2 其他求解常微分方程边值问题算法的简介 (6)5.3 模型求解 (8)5.3.1 运用MATLAB软件对模型求解结果及析 (8)6 模型的评价 (16)7 课程设计的总结与体会 (16)参考文献 (17)附录 (18)一、问题的提出1.1 问题背景：科学技术中常常需要求常微分方程的定解问题，微分方程里最简单的方程形式莫过于一阶常微分方程的初值问题，即：0(,)()dy f x y a x b dx y a y ⎧=≤≤⎪⎨⎪=⎩ （1）其中a ，b 为常数.虽然求解此类微分方程有各种各样的解析方法，但解析方法只能用于求解一些特殊类型方程，实际问题中归结出来的微分方程主要靠数值解法求解.因为一阶常微分方程简单但又是求解其他方程的基础，所以发展了许多典型的解法.本文着重讨论一类高精度的单步法——隆格库塔法，并且运用四阶的隆格库塔格式来求解初值问题，并且通过实例运用四阶的隆格库塔格式来求解初值问题，同时与显式与隐式的Euler 格式求解出的结果进行精度比较.1.2 问题的具体内容实例一：在区间[0,1]上采用经典的四阶隆格库塔方法求解微分方程1(0)1dy y x dx y ⎧=-++⎪⎨⎪=⎩，其精确解为x y x e -=+，步长为0.5,然后用欧拉法，预估校正法分别求解，且将计算结果与精确解进行比较，对三个算法的收敛性的进行分析比较.实例二：在区间[0,1]上用经典的四阶龙格库塔方法求解初值问题 (0)1x dy xe y dx y -⎧=-⎪⎨⎪=⎩， 其精确解为21(2)2xx e -+，然后用欧拉法，预估校正法分别求解，且将计算结果与精确解进行比较，对三个算法的收敛性的进行分析比较.最后在区间[0,1]上分别取步长h=0.1;0.05时进行计算，并且探究选取不同的步长对计算结果精度的影响.二、问题假设2.1 假设数值方法本身的计算是准确的.2.2 假设选取的步长趋于0时计算的结果会收敛到微分方程的准确解.2.3 假设步长的增加不会导致舍入误差的严重积累.三、符号系统3.1 符号说明符号含义h选取的步长*K平均斜率p精度的阶数∆前后两次计算结果的偏差y第n个节点的实验值n()y x第n个节点的精确值nδ实验值与精确值的绝对误差四、问题的分析问题要求运用隆格库塔算法来求解一阶微分方程的初值问题，针对前面提出的实例，本文先用经典的四阶隆格库塔法来求解上面的微分方程，为了体现隆格库塔法的优越性，同时用欧拉法，预估校正法分别求解，且将计算结果与精确解进行比较，对三个算法的收敛性的进行分析比较.最后在区间[0,1]上分别取步长h=0.1;0.05时进行计算，分析在选取不同的步长时，求解结果的精度变化如何.下面是欧拉法，预估校正法以及经典的四阶隆格库塔法的计算公式.4.1欧拉格式（1）显式欧拉格式1(,)n n n n y y hf x y +=+ （2） 局部截断误差22211()()()()22n n n h h y x y y y x o h ξ++''''-=≈= （3） （2）隐式欧拉格式111(,)n n n n y y hf x y +++=+ （4）局部截断误差2211()()()2n n n h y x y y x o h ++''-≈-= (5) 4.2 预估校正法预估 1(,)n n n n y y hf x y +=+ (6)校正 111[(,)(,)]2n n n n n n h y y f x y f x y +++=++ (7) 统一格式 1[(,)(,(,))]2n n n n n n n n h y y f x y f x h y hf x y +=++++ (8) 平均化格式 11(,),(,),1().2p n n n c n n p n p c y y hf x y y y hf x y y y y ++⎧⎪=+⎪=+⎨⎪⎪=+⎩ (9)4.3 四阶龙格库塔方法的格式（经典格式）112341213243(22),6(,),(,),22(,),22(,).n n n n n n n n n n 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 +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩(10)五、模型的建立与模型求解5.1 隆格库塔法的基本原理隆格库塔法是一种高精度的单步法，这类方法与下述Taylor 级数法有着紧密的联系.5.1.1 Taylor 级数设初值问题 '00(,)()y f x y y x y ⎧=⎨=⎩有解，按泰勒展开，有2'''1()()()()....2n n n n h y x y x hy x y x +=+++； （11） 其中()y x 的各阶导数依据所给方程可以用函数f 来表达，下面引进函数序列(,)j f x y 来描叙求导过程，即(0)(1)'(0)''(1)(1)(1)'''(2),f f y f f y f f x x f f y f f x y ∂∂=≡=+≡∂∂∂∂=+≡∂∂ （12）(2)(2)()(1)j j j j f f y f f x y ---∂∂=+≡∂∂ （13） 根据上式，结果导出下面 Taylor 格式2'''()1...2!!pp n n nn n h h y y hy y y p +=++++ （14）其中一阶Taylor 格式为： '1n n n y y hy +=+提高Taylor 格式的阶数p 即可提高计算结果的精度，显然，p 阶Taylor 格式的局部截断误差为：11(1)1(1)!p p n n h y y y p ζ+++-+=+ （15） 因此它具有p 阶精度.5.1.2 隆格库塔方法的基本思想隆格库塔法实质就是间接地使用Taylor 级数法的一种方法,考察差商1()()n n y x y x h+- 根据微分中值定理，这有'1()()()n n n y x y x y x h h θ+-=+ （16） 利用所给方程 '(,)y f x y =1()()(,())n n n n y x y x hf x h y x h θθ+=+++ （17） 设 平均斜率*(,())n n K f x h y x h θθ=++，由此可见，只要对平均斜率一种算法，便相应地可以导出一种计算格式.再考察改进的Euler 格式，它可以改写成平均化的形式：1121211()2(,)(,)n n n n n n h y y K K K f x y K f x y hK ++⎧=++⎪⎪=⎨⎪=+⎪⎩（18） 这个处理过程启示我们，如果设法在1(,)n n x x +内多预测几个点的斜率值，然后将它们加权平均作为平均斜率，则有可能构造具有更高精度的计算格式，这就是隆格库塔法的基本思想.5.1.3 四阶的隆格库塔法为了方便起见，本文主要运用经典的隆格库塔算法-四阶隆格库塔格式.其格式如下：112341213243(22),6(,),(,),22(,),22(,).n n n n n n n n n n 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 +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩（19）下面为其具体的算法流程图：5.2 其他求解常微分边值问题算法的简介5.2.1欧拉数值算法（显式）微分方程里最简单的方程形式莫过于一阶常微分方程的初值问题，即：(,)()dyf x y a x b dx y a y ⎧=≤≤⎪⎨⎪=⎩(20) 其中a ，b 为常数.开始输入x 0,y 0,h,N x 1=x 0+hk 1=f(x 0,y 0),k 2=f(x 0+h/2,y 0+hk 1/2)k 3=f(x 0+h/2,y 0+hk 2/2),k 4=f(x 1,y 0+hk 3)y 1=y 0+h(k 1+2k 2+2k 3+k 4)/6n=1输出x 1,y 1n=N? 结束n=n+1x 0=x 1,y 0=y 1否是图5.1 龙格-库塔法流程图因为其简单但又是求解其他方程的基础，所以发展了许多典型的解法.所有算法中的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 +=+ (21)简单欧拉法的算法过程介绍如下:给出自变量x 的定义域[,]a b ，初始值0y 及步长h .对0,1,()/k b a h =-，计算1(,)k k k k y y hf x y +=+5.2.2欧拉数值算法（隐式）隐式欧拉法也叫退欧拉法，隐式欧拉法的公式如下所示：111(,)n n n n y y hf x y +++=+ (22)隐式欧拉法是一阶精度的方法，比它精度高的公式是：111[(,)(,)]2n n n n n n hy y f x y f x y +++=++ (23)隐式欧拉的算法过程介绍如下.给出自变量x 的定义域[,]a b ，初始值0y 及步长h .对0,1,()/k b a h =-，用牛顿法或其他方法求解方111(,)k k k k y y hf x y +++=+得出1k y +.5.2.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 ++=+⎧⎪⎨=++⎪⎩ (24)四阶龙格-库塔法有多种形式，除了改进的欧拉法外还有中点法.中点法计算公式为：1[,(,)]22n n n n n n h h y y hf x y f x y +=+++ (25)5.3 模型求解5.3.1运用MATLAB 软件对模型求解结果及分析用欧拉法、改进的欧拉格式、经典的四阶龙格库塔法来求解常微分方程的边值问题，并且比较其精度（具体的MATLAB 源程序见附录） 以下进行实例分析：实例一. 1(0)1dyy x dx y ⎧=-++⎪⎨⎪=⎩由题可知精确解为：x y x e -=+，当x=0时，y(x)=0.在这里取步长h 为0.1， 通过MATLAB 程序的计算，相应的结果如下：表5-1 步长为0.1时各方法的预测值与精确值的比较(精确到6位小数)初值 Euler 法 相对误差 预估校正法 相对误差 经典四阶库 相对误差精确值 0 -- -- -- -- -- -- 1.00000 0.1 1.00910 0.00424 1.00500 0.00016 1.00484 0.00000 1.00484 0.2 1.02646 0.00759 1.01903 0.00029 1.01873 0.00000 1.01873 0.3 1.05134 0.01011 1.04122 0.00038 1.04082 0.00000 1.04082 0.4 1.08304 0.01189 1.07080 0.00045 1.07032 0.00000 1.07032 0.5 1.12095 0.01303 1.10708 0.00049 1.10653 0.00000 1.10653 0.6 1.16451 0.01366 1.14940 0.00052 1.14881 0.00000 1.14881 0.7 1.21319 0.01388 1.19721 0.00052 1.19659 0.00000 1.19659 0.8 1.26655 0.01378 1.24998 0.00052 1.24933 0.00000 1.24933 0.9 1.32414 0.01344 1.30723 0.00050 1.30657 0.00000 1.30657 1.01.38558 0.01294 1.36854 0.00048 1.36788 0.000001.36788步长为0.1时的精确值与预测值的比较精确值欧拉法改进欧拉格式四阶龙格库塔轴Y00.10.20.30.40.50.60.70.80.91 1.1 1.2 1.3 1.4X轴图5.2 步长为0.1时各方法的预测值与精确值的比较原函数图像轴YX轴图5.3 步长为0.1时原函数图像在这里取步长h为0.05，通过MATLAB程序的计算，相应的结果如下：表5-2 h=0.05时三个方法与精确值的真值表步长Euler法相对误差预估校正法相对误差经典四阶库相对误差精确值0 -- -- -- -- -- --1.00000 0.05 1.00250 0.00911 1.00125 0.01035 1.00123 0.01037 1.00484 0.10 1.00738 0.01711 1.00488 0.01954 1.00484 0.01958 1.01873 0.15 1.01451 0.02405 1.01076 0.02765 1.01071 0.02770 1.04082 0.20 1.02378 0.03001 1.01880 0.03473 1.01873 0.03479 1.07032 0.25 1.03509 0.03507 1.02889 0.04085 1.02880 0.04093 1.10653 0.30 1.04834 0.03930 1.04092 0.04610 1.04082 0.04619 1.14881 0.35 1.06342 0.04277 1.05480 0.05053 1.05469 0.05063 1.19659 0.40 1.08025 0.04555 1.07044 0.05422 1.07032 0.05432 1.24933 0.45 1.09874 0.04772 1.08775 0.05724 1.08763 0.05734 1.30657 0.50 1.11880 0.04933 1.10666 0.05964 1.10653 0.05975 1.36788 0.55 1.14036 0.05045 1.12709 0.06150 1.12695 0.06161 1.00484 0.60 1.16334 0.05113 1.14895 0.06286 1.14881 0.06298 1.01873 0.65 1.18768 0.05143 1.17219 0.06379 1.17205 0.06391 1.04082 0.70 1.21329 0.05139 1.19674 0.06433 1.19659 0.06445 1.07032 0.75 1.24013 0.05106 1.22252 0.06453 1.22237 0.06465 1.10653 0.80 1.26812 0.05048 1.24949 0.06443 1.24933 0.06455 1.14881 0.85 1.29722 0.04969 1.27757 0.06408 1.27742 0.06419 1.19659 0.90 1.32735 0.04871 1.30673 0.06349 1.30657 0.06361 1.249330.95 1.35849 0.04759 1.33690 0.06272 1.33674 0.06283 1.306571.00 1.39056 0.04634 1.36804 0.06178 1.36788 0.06189 1.36788欧拉格式 改进欧拉格式四阶龙格库塔精确值图5.4 步长为0.05时各方法的预测值与精确值的比较11.051.11.151.21.251.31.351.4X 轴Y 轴原函数图像图5.5 步长为0.05时原函数图像实例2 (0)1xdy xe ydx y -⎧=-⎪⎨⎪=⎩由题可知精确解为：21(2)2x x e -+当x=0时，y(x)=0.在这里取步长h为0.1，通过MATLAB程序的计算，相应的结果如下：步长Euler法相对误差预估校正法相对误差经典四阶库相对误差精确值0.1 0.9000 0.0318 0.9096 0.0214 0.9094 0.0216 0.9295 0.2 0.8192 0.0611 0.8359 0.0420 0.8356 0.0424 0.8726 0.3 0.7544 0.0876 0.7761 0.0614 0.7757 0.0619 0.8268 0.4 0.7027 0.1109 0.7277 0.0793 0.7272 0.0799 0.7903 0.5 0.6617 0.1310 0.6886 0.0956 0.6881 0.0963 0.7615 0.6 0.6294 0.1480 0.6572 0.1103 0.6567 0.1110 0.7387 0.7 0.6040 0.1621 0.6320 0.1234 0.6315 0.1241 0.7209 0.8 0.5841 0.1737 0.6116 0.1349 0.6111 0.1355 0.70690.9 0.5686 0.1830 0.5951 0.1450 0.5946 0.1456 0.69591.0 0.5563 0.1905 0.5815 0.1538 0.5811 0.1544 0.6872原函数图像轴Y00.10.20.30.40.50.60.70.80.91X轴图5.6 步长为0.1时原函数图像各方法的预测值与精确值的比较欧拉格式改进的格式四阶龙格库塔精确值轴Y0.10.20.30.40.50.60.70.80.91X轴图5.7 步长为0.1时各方法的预测值与精确值的比较在这里取步长h为0.05，通过MATLAB程序的计算，相应的结果如下：表5-4 步长为0.1时各方法的预测值与精确值的比较(精确到5位小数)步长Euler法相对误差预估校正法相对误差经典四阶库相对误差精确值0.050.95000 0.01342 0.95245 0.01088 0.95242 0.01090 0.96292 0.100.90490 0.02650 0.90947 0.02158 0.90942 0.02163 0.92953 0.150.86428 0.03916 0.87067 0.03206 0.87061 0.03213 0.89951 0.200.82774 0.05137 0.83567 0.04228 0.83559 0.04237 0.87256 0.250.79491 0.06307 0.80414 0.05219 0.80405 0.05230 0.84842 0.300.76545 0.07423 0.77576 0.06176 0.77566 0.06188 0.82682 0.350.73904 0.08482 0.75023 0.07096 0.75013 0.07110 0.80754 0.40 0.71541 0.09482 0.72730 0.07977 0.72719 0.07991 0.79035 0.45 0.69427 0.10422 0.70672 0.08817 0.70660 0.08832 0.77505 0.50 0.67539 0.11302 0.68825 0.09615 0.68813 0.09630 0.76146 0.55 0.65855 0.12123 0.67168 0.10370 0.67156 0.10386 0.74940 0.60 0.64352 0.12886 0.65683 0.11084 0.65671 0.11100 0.73871 0.65 0.63012 0.13593 0.64351 0.11756 0.64340 0.11773 0.72925 0.70 0.61819 0.14245 0.63157 0.12388 0.63145 0.12405 0.72087 0.75 0.60754 0.14847 0.62085 0.12981 0.62073 0.12998 0.71347 0.80 0.59805 0.15400 0.61121 0.13537 0.61110 0.13553 0.70691 0.85 0.58957 0.15908 0.60254 0.14058 0.60243 0.14073 0.70109 0.90 0.58198 0.16374 0.59470 0.14545 0.59460 0.14560 0.695930.95 0.57517 0.16802 0.58761 0.15001 0.58751 0.15016 0.691321.00 0.56904 0.17194 0.58117 0.15429 0.58107 0.15443 0.687190.10.20.30.40.50.60.70.80.910.550.60.650.70.750.80.850.90.951X 轴Y 轴各方法下的预测值与精确值的比较欧拉格式改进欧拉格式四阶龙格库塔精确值图5.8 步长为0.05时各方法的预测值与精确值的比较六、模型的评价本文着重讨论了4阶的隆格库塔法来求解微分方程，并且通过两个实例验证了隆格库塔法在求解初值问题的优越性.从上面的实例比较可知，在计算精度上，四阶经典龙格-库塔方法的误差最小，改进欧拉方法其次，欧拉方法误差则比较大，所以四阶经典龙格-库塔方法得到最佳的精度.而在计算量上面，相应地，很明显的四阶经典龙格-库塔方法也是最大，改进欧拉方法其次，欧拉方法计算量最小.这样的结果，说明了运用以上三种方法时，其计算量的多少与精度的大小成正比.我们在实际运用与操作中，可以根据实际情况，选择这3种方法中的其中一种最适合的，追求精度的话，可以使用四阶经典龙格-库塔方法；而改进的欧拉方法，在精度上和计算量上都表现得很出色，能够满足一般情况；而欧拉方法更主要的是适用于对y的估计上，而精度则有所欠缺，以上各方法的选择，都取决于具体的情况.七、课程设计的总结与体会本文着重采用隆格库塔法运用MATLAB编程来求解微分方程，相比于欧拉法以及预估校正法，隆格库塔法在提高近似值解的精度上是非常起作用的，而且又具有计算量不大、算法组织容易.其次，每一次的课程设计总是让我学到了更多的知识，不论是C++、SPASS还是MATLAB软件，这些都让我学到了如何解决实际问题的好工具，通过这些工具，是自己能够得到突破和成长.以下是我完成此次课程设计的几点体会：（1）必须学好基础知识，在做的过程中，发现自己有很多东西都不懂，要博学必须从一点一点做起.以往训练得少只是把握的不牢靠.所以做起来感到有点吃力.所以，无论什么学科，一定要打好基础.（2）程序设计要靠多练，多见识，那样形成一种编程思维，我想对我是有很大好处的.尤其像我这种平时学得不扎实的人.（3）做事情要有恒心，遇到困难不要怕，坚决去做.如果做出来了，固然高兴，如果没有做出来也没关系，自己努力了对得起自己就好.同时，把它看做是对自己的锻炼. （4）做程序特别是做大程序是很有趣的.虽然有的问题很难，要花很多时间很多精力，但是那种解决了一个问题时的喜悦足以把付出的辛苦补偿回来.得到一种心里的慰藉.参考文献[1] 李庆杨，王能超，易大义编．数值分析（第四版）[M]：华中科技大学出版社，2006.[2] 姜启源，谢金星，叶俊编．数学模型（第三版）[M]．北京：高等教育出版社，2005[3] 刘琼荪，数学实验[M]，高等教育出版社，2004[4] 王建伟，MATLAB7.X程序设计[M]，中国水利水电出版社，2007[5] 王高雄，周之铭等编.常微分方程（第三版） [M]:高等教育出版社，2006[6]何坚勇编著. 运筹学基础（第二版）[M]. 北京：清华大学出版社，2008.附录附录一：显示欧拉法matlab程序%欧拉法clear allclcx=[];y=[];y1=[];h=0.1;x=0:h:1;n=length(x);for i=1:ny(i)=f1(x(i));endfigure(1)plot(x,y,'g-');hold ony1(1)=1;for j=2:ny1(j)=y1(j-1)+h*f(x(j-1),y1(j-1)); endY=[x;y1];fid=fopen('data.txt','wt');fprintf(fid,'%6.2f %12.4f\n',Y);fclose(fid);plot(x,y1,'r-');figure(2)DT=abs(y-y1);plot(x,DT)%1.建立导数函数文件function z=f(x,y)z=y-2*x/y;%2.建立原函数文件function z1=f1(x)z1=(2*x+1)^(1/2);迭代n次的后退的欧拉格式matlab程序%1.建立导数函数文件function z=f(x,y)z=y-2*x/y;%2.建立原函数文件function z1=f1(x)z1=(2*x+1)^(1/2);%迭代n次的后退的欧拉格式clear allclcN=input('请输入迭代次数:');x=[];y=[];y1=[];h=0.1;x=0:0.1:1;n=length(x);for i=1:ny(i)=f1(x(i));endfigure(1)plot(x,y,'g-');for j=2:ny1(j)=y1(j-1)+h*f(x(j-1),y1(j-1)); T=y1(j);for k=1:NT=y1(j-1)+h*f(x(j),T);endy1(j)=T;endY=[x;y1];fid=fopen('data.txt','wt');fprintf(fid,'%6.2f %12.4f\n',Y); fclose(fid);plot(x,y1,'r-');figure(2)DT=abs(y-y1);plot(x,DT)附录二：预估校正法matlab程序%1.建立导数函数文件function z=f(x,y)z=y-2*x/y;%2.建立原函数文件function z1=f1(x)z1=(2*x+1)^(1/2);%预估校正法%欧拉法clear allclcx=[];y1(1)=1;y2=[];y2(1)=1;h=0.1;x=0:0.1:1;n=length(x);for i=1:ny(i)=f1(x(i));endfigure(1)plot(x,y,'g-');hold onfor j=2:ny1(j)=y2(j-1)+h*f(x(j-1),y2(j-1));y2(j)=y2(j-1)+(h/2).*[f(x(j-1),y2(j-1))+f(x(j),y1(j))]; endY=[x;y2];fid=fopen('data.txt','wt');fprintf(fid,'%6.2f %12.4f\n',Y);fclose(fid);plot(x,y2,'r-');figure(2)DT=y-y2;plot(x,DT)附录三：四阶龙格库塔法matlab程序%四阶龙格库塔法clear allclcn=length(x);y=[];y1=[];y1(1)=1;for i=1:ny(i)=f1(x(i));endfigure(1)plot(x,y,'g-');hold onfor j=2:nK1=f(x(j-1),y1(j-1));K2=f(x(j-1)+h/2,y1(j-1)+h/2*K1);K3=f(x(j-1)+h/2,y1(j-1)+h/2*K2);K4=f(x(j-1)+h,y1(j-1)+h*K3);y1(j)=y1(j-1)+(h/6)*(K1+2*K2+2*K3+K4); endY=[x;y1];fid=fopen('data1.txt','wt');fprintf(fid,'%6.2f %12.4f\n',Y);fclose(fid);plot(x,y1,'r-');figure(2)DT=abs(y-y1);plot(x,DT)。
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版
微分方程组的龙格-库塔公式求解matlab 版南京大学王寻1.一阶常微分方程组考虑方程组()()()()⎩⎨⎧====0000z x z ,z ,y ,x g 'z y x y ,z ,y ,x f 'y 其经典四阶龙格-库塔格式如下：对于n =0,1,2,...,计算()()⎪⎩⎪⎨⎧++++=++++=++4321143211226226L L L L h z z K K K K h y y n n n n 其中()()()()⎪⎪⎪⎩⎪⎪⎪⎨⎧+++=+++=⎪⎭⎫ ⎝⎛+++=⎪⎭⎫ ⎝⎛+++=⎪⎭⎫ ⎝⎛+++=⎪⎭⎫ ⎝⎛+++===33433422322311211211222222222222hL z ,hK y ,h x g L ,hL z ,hK y ,h x f K hL z ,hK y ,h x g L ,hL z ,hK y ,h x f K hL z ,hK y ,h x g L ,hL z ,hK y ,h x f K z ,y ,x g L ,z ,y ,x f K n n n n n n n n n n n n n n n n n n n n n n n n 下面给出经典四阶龙格-库塔格式解常微分方程组的matlab 通用程序：%marunge4s.m%用途：4阶经典龙格库塔格式解常微分方程组y'=f(x,y),y(x0)=y0%格式：[x,y]=marunge4s(dyfun,xspan,y0,h)%dyfun 为向量函数f(x,y),xspan 为求解区间[x0,xn],%y0为初值向量，h 为步长，x 返回节点，y 返回数值解向量function [x,y]=marunge4s(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y=zeros(length(y0),length(x));y(:,1)=y0(:);for n=1:(length(x)-1)k1=feval(dyfun,x(n),y(:,n));k2=feval(dyfun,x(n)+h/2,y(:,n)+h/2*k1);k3=feval(dyfun,x(n)+h/2,y(:,n)+h/2*k2);k4=feval(dyfun,x(n+1),y(:,n)+h*k3);y(:,n+1)=y(:,n)+(h/6).*(k1+2*k2+3*k3+k4);end如下为例题：例1：取h =0.02，利用程序marunge4s.m 求刚性微分方程组()()⎩⎨⎧=-==--=10100209999010z ,z 'z ,y ,z .y .'y 的数值解，其解析解为：xx x .e z ,e e y 100100010---=+=解：首先编写M 函数dyfun.m%dyfun.m function f=dyfun(t,y)f(1)=-0.01*y(1)-99.99*y(2);f(2)=-100*y(2);f=f(:);然后编写一个执行函数：close all ;clear all ;clc;[x,y]=marunge4s(@dyfun,[0500],[21],0.002);plot(x,y);axis([-50500-0.52]);text(120,0.4,'y(x)');text(70,0.1,'z(x)');title('数值解');figure;y1=exp(-0.01*x)+exp(-100*x);%解析解，用来做对比的。
Matlab在常微分方程求解中的应用
求解常微分方程初值问题的数值解的整体思路：
dy f (x, y) dx y(x0 ) y0
（1）
寻求准确解 y(x ) 在一系列离散节点：
x0 x1 x2 xn 上的近似值 y0 , y1, y2 ,, yn ,
yn 称为问题的数值解，数值解所满足的离散方程统称为差
分格式， hi xi xi1 称为步长，实用中常取定步长。
ode23t ode23tb
Trapezoidal rule with a “free” interpolant
Implicit Runge-Kutta formula with a backward differentiation formula of order two
Stiff ODE Stiff ODE Stiff ODE
2、使用Matlab软件求数值解时，高阶微分方程必须 等价地变换成一阶微分方程组.
Van der Pol微分方程
1928年荷兰的范德波耳(Van der Pol)为描述LC 回路的电子管振荡器建立了著名的vanderPol方 程.它在自激振荡理论中有着重要的意义,一直作 为数学物理方程中的一个基本方程.这是一个具有 可变非线性阻尼的微分方程,代表了一类极为典型 的非线性问题.和其他非线性微分方程在数学上无 法精确求解一样,人们一直在努力寻找求解这类方 程近似解析解的方法,并乐于用Van der Pol方程 来检验求解方法的有效性.
>> type vdp1.m function dy = vdp1(t, y) mu = 1; dy = [y(2); mu*(1-y(1)^2)*y(2)-y(1)];
建立了vdp1.m后，即可选用前述ODE指令来求解
在 =1时，van der Pol方程并非 Stiff系统，所 以使用ode45来显示结果 odeBasic01.m
matlab用四阶龙格库塔函数求解微分方程组
一、介绍Matlab作为一种强大的科学计算软件，提供了众多函数和工具来解决微分方程组。
其中，四阶龙格库塔函数是一种常用的数值方法，用于求解常微分方程组。
本文将介绍如何使用Matlab中的四阶龙格库塔函数来求解微分方程组，并对该方法的原理和实现进行详细说明。
二、四阶龙格库塔方法四阶龙格库塔方法是一种常用的数值方法，用于求解常微分方程组。
它是一种显式的Runge-Kutta方法，通过逐步逼近微分方程的解，在每一步使用多个中间值来计算下一步的解。
该方法通过四个中间值来计算下一步的状态，并且具有较高的精度和稳定性。
三、在Matlab中使用四阶龙格库塔方法求解微分方程组在Matlab中，可以使用ode45函数来调用四阶龙格库塔方法来解决微分方程组的问题。
ode45函数是Matlab提供的用于求解常微分方程组的函数，可以通过指定微分方程组以及初值条件来调用四阶龙格库塔方法来进行求解。
1. 定义微分方程组我们需要定义要求解的微分方程组。
可以使用Matlab中的匿名函数来定义微分方程组，例如：```matlabf = (t, y) [y(2); -sin(y(1))];```其中，f是一个匿名函数，用于表示微分方程组。
在这个例子中，微分方程组是y' = y2, y2' = -sin(y1)。
2. 指定初值条件和求解区间接下来，我们需要指定微分方程组的初值条件和求解区间。
初值条件可以通过指定一个初始时刻的状态向量来完成，例如：```matlabtspan = [0, 10];y0 = [0, 1];```其中，tspan表示求解区间，y0表示初值条件。
3. 调用ode45函数进行求解我们可以通过调用ode45函数来求解微分方程组的数值解。
具体的调用方式如下：```matlab[t, y] = ode45(f, tspan, y0);```其中，t和y分别表示求解的时间点和对应的状态值。
四、示例下面我们通过一个具体的例子来演示如何使用Matlab中的四阶龙格库塔方法来求解微分方程组。
matlab编的4阶龙格库塔法解微分方程的程序
matlab编的4阶龙格库塔法解微分方程的程序2010-03—10 20：16function varargout=saxplaxliu(varargin)clc，clearx0=0；xn=1。
2；y0=1；h=0。
1;[y,x］=lgkt4j（x0,xn，y0,h)；n=length（x)；fprintf（’ i x（i) y(i)\n’)；for i=1:nfprintf('%2d %4.4f ％4.4f\n’,i，x(i），y（i))；endfunction z=f（x，y）z=-2*x*y^2；function [y,x]=lgkt4j（x0，xn,y0，h)x=x0:h:xn；n=length（x）；y1=x;y1（1)=y0;for i=1:n-1K1=f(x(i)，y1(i））;K2=f（x（i)+h/2,y1（i)+h/2*K1)；K3=f(x(i)+h/2,y1（i)+h/2*K2);K4=f(x（i）+h,y1（i）+h*K3）;y1(i+1)=y1(i）+h/6＊(K1+2*K2+2＊K3+K4）;endy=y1；结果：i x（i） y（i)1 0.0000 1。
00002 0。
1000 0.99013 0。
2000 0.96154 0.3000 0。
91745 0.4000 0。
86216 0.5000 0。
80007 0。
6000 0。
73538 0。
7000 0.67119 0。
8000 0。
609810 0。
9000 0。
552511 1.0000 0.500012 1.1000 0.452513 1.2000 0.4098yi+1=yi+h＊（ K1+ K2)/2K1=f（xi,yi)K2=f(xi+h，yi+h*K1)依次类推，如果在区间［xi，xi+1］内多预估几个点上的斜率值K1、K2、……Km,并用他们的加权平均数作为平均斜率K*的近似值，显然能构造出具有很高精度的高阶计算公式。
四阶Runge-Kutta法求解常微分方程
【SOS】在matlab中四阶Runge-Kutta法求解常微分方程匚悬赏分：150 -解决时间：2008-9-16 10:19【SOS】在matlab中四阶Runge-Kutta法求解常微分方程dx1/dt=x2 x1(0)=1e-8dx2/dt=x3 x2(0)=0dx3/dt=x4 x3(0)=0dx4/dt=f(t,x1,x2,x3,x4) x4(0)=0求解区间［0,1e-6］，在matlab中用四阶Runge-Kutta法求解，怎么用从一阶循环到四阶？有没有类似的具体算例？最好能帮小弟写个子程序希望各路高人出手相助，关乎小弟能否明年一月份顺利毕业，在此先多谢各位高人了。
提问者：zhihaikou -三级最佳答案没试过matlab，算这玩意太慢了，有fortran版的要不，有兴趣的话可以参考一下。
代码:SUBROUTINE run ge_kutta()!关于Runge-Kutta方法，该方法是用来解形如y'=f(t,y)的常微分方程的!经典的4阶R-K方法的公式如下：! Yn+1 = Yn + h/6 * (K1+2K2+2K3+K4)!其中! K仁f(Tn,Yn)! K2=f(T n+h/2,Y n+h/2*K1)! K3=f(T n+h/2,Y n+h/2*K2)! K4=f(T n+h,Y n+h*K3)type( no de_fish),po in ter :: pi7HO!S%（『!）siuo!」inu% 厂u 人 NOa%（r!）siuau;nu%r u 人NOd%（『!）siuo!」inu% 厂 u 人t7HN%（『!）siuo!」inu% 厂 u 人£ON%（『!）siuo!」inu%厂u 人（（>HO!S%lueu;nu _U!UJ >（ NOa%lueu;nu _u!iu > （ NOd%lueu;nu _u!iu >（ >HN%lueu;nu _U!UJ > eON%lueu;nu _u!iu >WW : >d31S -一ii（W+£>k0N+Z>k0N+"）*（09dois —|j ） + u 人二厂u 人L+u 人直具 *-以吐u 人甲：ed3丄S-—一ii（e>i*de;s -p+u 人'dois —：M+iuaurK ）—：|j ）uo!Qun 厂 |E 」60屮!=厂》£>UL|+u 人 POM u人'u+l DOMfeJB來：kSd3丄S ------- i（2>1*（0S/deis _p ）+UA 10S/deis _p+juajjno _p ）uo!pun4_|ej6aju!=e>j◊u ⑵q ）+u 人［1 變歪 u 人'乙/i|+i pOMfeJB Ux 來：c sd3丄s ------------- i（L>1*（0S/deis _|j ）+UA 107/deis _p+juajjno _p ）uo!pun4_|ej6aju!=3>j以*⑵L|）+u 人 H 變歪U 人'乙/1|+1 POMHB 辽X 來：2 2 C □丄s --------- i（uA^uejjno-pJuoipunPie 」6。
matlab算法-求解微分方程数值解和解析解
MATLAB是一种用于数学计算、工程和科学应用程序开发的高级技术计算语言和交互式环境。
它被广泛应用于各种领域，尤其在工程和科学领域中被用于解决复杂的数学问题。
微分方程是许多工程和科学问题的基本数学描述，求解微分方程的数值解和解析解是MATLAB算法的一个重要应用。
1. 求解微分方程数值解在MATLAB中，可以使用各种数值方法来求解微分方程的数值解。
其中，常见的方法包括欧拉法、改进的欧拉法、四阶龙格-库塔法等。
这些数值方法可以通过编写MATLAB脚本来实现，从而得到微分方程的近似数值解。
以常微分方程为例，可以使用ode45函数来求解微分方程的数值解。
该函数是MATLAB中用于求解常微分方程初值问题的快速、鲁棒的数值方法，可以有效地得到微分方程的数值解。
2. 求解微分方程解析解除了求解微分方程的数值解外，MATLAB还可以用于求解微分方程的解析解。
对于一些特定类型的微分方程，可以使用符号计算工具箱中的函数来求解微分方程的解析解。
通过符号计算工具箱，可以对微分方程进行符号化处理，从而得到微分方程的解析解。
这对于研究微分方程的性质和特点非常有帮助，也有助于理论分析和验证数值解的准确性。
3. MATLAB算法应用举例在实际工程和科学应用中，MATLAB算法求解微分方程问题非常常见。
在控制系统设计中，经常需要对系统的动态特性进行分析和设计，这通常涉及到微分方程的建模和求解。
通过MATLAB算法，可以对系统的微分方程进行数值求解，从而得到系统的响应曲线和动态特性。
另外，在物理学、生物学、经济学等领域的建模和仿真中，也经常需要用到MATLAB算法来求解微分方程问题。
4. MATLAB算法优势相比于其他数学软件和编程语言，MATLAB在求解微分方程问题上具有明显的优势。
MATLAB提供了丰富的数值方法和工具，能够方便地对各种微分方程进行数值求解。
MATLAB具有直观的交互式界面和强大的绘图功能，能够直观地展示微分方程的数值解和解析解，有利于分析和理解问题。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
姓名：樊元君学号：2012200902 日期：2012.11.06 一、实验目的
掌握MATLAB语言、C/C++语言编写计算程序的方法、掌握改进欧拉法与四阶龙格-库塔求解一阶常微分方程的初值问题。
掌握使用MATLAB程序求解常微分方程问题的方法。
二、实验内容
1、分别写出改进欧拉法与四阶龙格-库塔求解的算法，编写程序上机调试出结果，要求所编程序适用于任何一阶常微分方程的数值解问题，即能解决这一类问题，而不是某一个问题。
实验中以下列数据验证程序的正确性。
求，步长h=0.25。
2、实验注意事项
的精确解为，通过调整步长，观察结果的精度的变化
三、程序流程图：
●改进欧拉格式流程图：
●四阶龙格库塔流程图：
四、源程序：
●改进后欧拉格式程序源代码：
format long
h=input('h=');
x0=input('x0=');
y0=input('y0=');
disp('输入的范围是：');
X=input('X=');Y=input('Y=');
n=round((Y-X)/h);
i=1;x1=0;yp=0;yc=0;
for i=1:1:n
x1=x0+h;
yp=y0+h*(-x0*(y0)^2);%yp=y0+h*(y0-2*x0/y0);% yc=y0+h*(-x1*(yp)^2);%yc=y0+h*(yp-2*x1/yp);% y1=(yp+yc)/2;
x0=x1;y0=y1;
y=2/(1+x0^2);%y=sqrt(1+2*x0);%
fprintf('结果=%.3f,%.8f,%.8f\n',x1,y1,y);
end
end
●四阶龙格库塔程序源代码：
format long
h=input('h=');
x0=input('x0=');
y0=input('y0=');
disp('输入的范围是：');
X=input('X=');Y=input('Y=');
n=round((Y-X)/h);
i=1;x1=0;k1=0;k2=0;k3=0;k4=0;
for i=1:1:n
x1=x0+h;
k1=-x0*y0^2;%k1=y0-2*x0/y0;%
k2=(-(x0+h/2)*(y0+h/2*k1)^2);%k2=(y0+h/2*k1)-2*(x0+h/2)/(y0+h/2*k1);% k3=(-(x0+h/2)*(y0+h/2*k2)^2);%k3=(y0+h/2*k2)-2*(x0+h/2)/(y0+h/2*k2);% k4=(-(x1)*(y0+h*k3)^2);%k4=(y0+h*k3)-2*(x1)/(y0+h*k3);%
y1=y0+h/6*(k1+2*k2+2*k3+k4);%y1=y0+h/6*(k1+2*k2+2*k3+k4);%
x0=x1;y0=y1;
y=2/(1+x0^2);%y=sqrt(1+2*x0);%
fprintf('结果=%.3f,%.7f,%.7f\n',x1,y1,y);
end
end
五、运行结果：
改进欧拉格式结果：
四阶龙格库塔结果：
步长分别为：0.25和0.1时，不同结果显示验证了步长减少，对于精度的提高起到很大作用，有效数字位数明显增加。
六、实验小结：
通过这次实验学习，首先第一点对改进欧拉格式和四阶龙格库塔的原理推导有了深入的理解，改进欧拉格式采用（预报+校正）模式得到较精确的原函数数值解；而四阶龙格库塔则采用多预报几个点的斜率值，采用加权平均作为平均斜率的近似值的思想达到更高精度的数值解，二阶龙格库塔的特例就是改进后的欧拉格式。
七、思考题：
如何对四阶龙格-库塔法进行改进，以保证结果的精度。
答：可以通过计算结果的精度处理步长来保证结果的精度。
（1）步长折半。
对于给定精度ε，如果某次计算结果精度(/2)11||h h i i y y ++∆=-，ε∆>，反复将步长减半，直到ε∆<，这时的(/2)1h i y +可作为结果。
（2）步长加倍。
对于给定精度ε，如果某次计算结果精度(/2)11||h h i i y y ++∆=-，ε∆<，反复将步长加倍，直到ε∆>，这时的()1h i y +可作为结果。
Welcome To Download !!!
欢迎您的下载，资料仅供参考！。