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*五、运行结果:改进欧拉格式结果:;}四阶龙格库塔结果:步长分别为:和时,不同结果显示验证了步长减少,对于精度的提高起到很大作用,有效数字位数明显增加。
matlab龙格库塔方法求解二元二阶常微分方程组

matlab龙格库塔方法求解二元二阶常微分方程组文章标题:深入探讨matlab中的龙格库塔方法及其在求解二元二阶常微分方程组中的应用摘要:在科学与工程领域,常常需要求解复杂的微分方程组,而matlab作为一种强大的数学工具,提供了许多求解微分方程组的方法。
本文将深入探讨matlab中的龙格库塔方法及其在求解二元二阶常微分方程组中的应用,以便读者全面理解该方法并能灵活应用于实际问题中。
正文:一、介绍龙格库塔方法龙格-库塔法(Runge-Kutta methods)是一种数值求解常微分方程的方法,通过将微分方程的解进行离散化,将微分方程转化为差分方程,从而进行数值求解。
龙格库塔方法通过迭代计算,能够得到微分方程的数值解,广泛应用于科学计算和工程技术领域。
二、matlab中的龙格库塔方法在matlab中,龙格库塔方法通过ode45函数实现,该函数能够对一阶或高阶常微分方程进行数值求解。
用户可以通过设定初始条件、微分方程表达式,以及积分区间等参数,快速得到微分方程的数值解。
ode45函数采用自适应步长的方式进行求解,能够有效解决微分方程解的数值稳定性和精确度问题。
三、龙格库塔方法在求解二元二阶常微分方程组中的应用考虑如下形式的二元二阶常微分方程组:$$\begin{cases}y_1' = f_1(t, y_1, y_2) \\y_2' = f_2(t, y_1, y_2)\end{cases}$$其中,$y_1(t)$和$y_2(t)$是未知函数,$f_1(t, y_1, y_2)$和$f_2(t,y_1, y_2)$分别表示其对应的函数表达式。
通过matlab中的ode45函数,可以将该二元二阶常微分方程组转化为一阶常微分方程组的形式,然后利用龙格库塔方法进行数值求解。
设定初始条件$y_1(0) = y1_0, y_2(0) = y2_0$,对应的一阶方程组为:$$\begin{cases}u_1' = u_3 \\u_2' = u_4 \\u_3' = f_1(t, u_1, u_2) \\u_4' = f_2(t, u_1, u_2)\end{cases}$$其中,$u_1(t) = y_1(t), u_2(t) = y_2(t), u_3(t) = y_1'(t), u_4(t) =y_2'(t)$,通过ode45函数求解该一阶常微分方程组即可得到原二元二阶常微分方程组的数值解。
微分方程的数值解法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常微分⽅程数值解——欧拉法、改进的欧拉法与四阶龙格库塔⽅法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。
matlab求解常微分方程组

matlab求解常微分方程组常微分方程组是数学中的一个重要分支,它描述了多个变量随时间变化的关系。
在实际应用中,常微分方程组经常被用来描述物理、化学、生物等领域中的动态系统。
本文将介绍如何使用MATLAB求解常微分方程组。
MATLAB是一种强大的数学软件,它提供了许多工具和函数来求解常微分方程组。
在MATLAB中,我们可以使用ode45函数来求解常微分方程组。
ode45函数是一种常用的数值求解器,它使用龙格-库塔方法来求解常微分方程组。
我们需要定义常微分方程组。
常微分方程组通常采用向量形式表示,例如:dy/dt = f(t,y)其中,y是一个向量,f(t,y)是一个向量函数。
在MATLAB中,我们可以使用匿名函数来定义f(t,y)。
例如,如果我们要求解以下常微分方程组:dy1/dt = -y1 + 2*y2dy2/dt = -2*y1 + 3*y2我们可以定义f(t,y)为:f = @(t,y) [-y(1) + 2*y(2); -2*y(1) + 3*y(2)];接下来,我们需要指定初值条件。
初值条件是指在t=0时,y的值。
在MATLAB中,我们可以使用一个向量来表示初值条件。
例如,如果我们要求解以下常微分方程组:dy1/dt = -y1 + 2*y2dy2/dt = -2*y1 + 3*y2初值条件为:y(0) = [1; 0]我们可以定义初值条件为:y0 = [1; 0];现在,我们可以使用ode45函数来求解常微分方程组。
ode45函数的语法如下:[t,y] = ode45(f,tspan,y0)其中,f是一个函数句柄,tspan是一个包含起始时间和结束时间的向量,y0是一个包含初值条件的向量。
ode45函数将返回一个包含时间和解向量的矩阵。
例如,如果我们要求解以下常微分方程组:dy1/dt = -y1 + 2*y2dy2/dt = -2*y1 + 3*y2初值条件为:y(0) = [1; 0]时间范围为0到10秒,我们可以使用以下代码来求解:f = @(t,y) [-y(1) + 2*y(2); -2*y(1) + 3*y(2)];tspan = [0 10];y0 = [1; 0];[t,y] = ode45(f,tspan,y0);现在,我们可以绘制解向量随时间变化的图像。
matlab四阶龙格库塔法解方程组

matlab四阶龙格库塔法解方程组摘要:一、引言二、龙格库塔法介绍1.龙格库塔法的基本原理2.龙格库塔法的发展历程三、MATLAB 实现四阶龙格库塔法1.MATLAB 中龙格库塔法的函数2.四阶龙格库塔法的MATLAB 实现四、龙格库塔法解方程组的应用1.线性方程组的求解2.非线性方程组的求解五、结论正文:一、引言在数学领域,求解方程组是一项基本任务。
龙格库塔法作为高效数值求解线性方程组的方法,被广泛应用于各个领域。
MATLAB 作为一款强大的数学软件,可以方便地实现龙格库塔法求解方程组。
本文将介绍MATLAB 中四阶龙格库塔法解方程组的原理、实现与应用。
二、龙格库塔法介绍1.龙格库塔法的基本原理龙格库塔法(Runge-Kutta method)是一种求解常微分方程初值问题的数值方法。
它通过求解一组线性方程来逼近微分方程的解,具有较高的数值稳定性和精度。
龙格库塔法可以分为四阶、五阶等多种形式,其中四阶龙格库塔法是较为常用的一种。
2.龙格库塔法的发展历程龙格库塔法由德国数学家卡尔·龙格(Carl Runge)和英国数学家詹姆斯·库塔(James Kutta)分别在1900 年和1901 年独立发现。
自那时以来,龙格库塔法在数学、物理、工程等领域得到了广泛应用,并发展出了多种改进和扩展。
三、MATLAB 实现四阶龙格库塔法1.MATLAB 中龙格库塔法的函数在MATLAB 中,可以使用内置函数ode45、ode23、ode113 等实现龙格库塔法求解常微分方程。
这些函数分别对应四阶、五阶和三阶龙格库塔法。
2.四阶龙格库塔法的MATLAB 实现以下是一个使用MATLAB 实现四阶龙格库塔法求解方程组的示例:```matlabfunction [x, status] = solve_system_with_ode45(A, B, x0)% 定义方程组func = @(t, x) A * x + B;% 初始条件x0 = [1; 2];% 时间区间tspan = [0, 10];% 求解[x, status] = ode45(func, tspan, x0);end```四、龙格库塔法解方程组的应用1.线性方程组的求解线性方程组在数学、物理、工程等领域具有广泛应用。
matlab迭龙格库塔法解常微分方程

一、介绍迭龙格-库塔法(Runge-Kutta method)是一种数值求解常微分方程(ODE)的常用方法。
它是由卡尔·迭龙格(Carl Runge)和马丁·威尔黑尔姆·库塔(Wilhelm Kutta)在20世纪初提出的,该方法以两位数值分析家的名字来命名。
二、简单描述迭龙格-库塔法是通过数值逼近的方式,来计算常微分方程的近似解。
它是一种显式求解方法,适用于解非线性常微分方程和具有较大阶数的常微分方程。
三、数学原理迭龙格-库塔法主要是通过将微分方程转化为差分方程,利用数值解的方式来逼近微分方程的解。
它是一种显式方法,通过不断迭代得到下一个时间步的近似解。
四、matlab中的应用在matlab中,可以使用ode45函数来调用迭龙格-库塔法求解常微分方程。
ode45函数是matlab中集成的一个函数,通过调用ode45函数,可以直接求解常微分方程的数值解。
五、实例演示下面通过一个简单的例子来演示如何使用matlab中的ode45函数来求解常微分方程。
我们考虑一个简单的一阶常微分方程:dy/dt = -y初始条件为y(0) = 1。
在matlab中,可以通过以下代码来求解该微分方程:```定义微分方程的函数function dydt = myode(t, y)dydt = -y;调用ode45函数求解[t, y] = ode45(myode, [0, 5], 1);plot(t, y);```运行以上代码,即可得到微分方程的数值解,并通过绘图来展示解的变化。
六、总结迭龙格-库塔法是一种常用的数值解常微分方程的方法,它在matlab中有较为方便的调用方式。
通过ode45函数,可以快速求解常微分方程的数值解,并通过绘图来展示结果。
希望本篇文章对读者有所帮助,谢谢阅读。
七、应用场景和优势在实际应用中,迭龙格-库塔法广泛应用于各种科学和工程领域,如物理学、化学、生物学、经济学等。
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微分方程的解法

-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
Matlab求解教程:求解延迟微分方程

Matlab求解教程:求解延迟微分方程使用MA TLAB求解延时微分方程的两种方法:DDE23和SimuLink有些不同点需要注意,否则结果会出现错误使用MA TLAB来求解延迟微分方程是在生物数学和化学计算求解中经常遇到的事,在其它领域也比较常见。
我所知道的,在MA TLAB中求解微分方程有三种方法:1.使用ode45(龙格-库塔法的一个变种)求解,这时用一个数组,记录y的延迟项,但是c的值要考虑步长,再代入方程就能实现延时效应;2.使用dde23求解常数延时方程、使用ddesd可以用来求解延迟与时间t有关的延迟微分方程;3.使用SimuLink建模求解,SimuLink是求解广义微分代数系统的通用工具,功能很强大,但是看惯了编程指令的人可能不大习惯,调试似乎也不太方便。
既然本文专门讨论求解延迟微分方程,就先介绍一下专用工具dde23,dde系列求解函数是由Southern Methodist University 的L.F. Shampine 和S. Thompson根据他们早期使用Fortran编写的Fortran 90 DDE Solver 移植到MATLAB上的,从MA TLAB6.5开始加入MA TLAB的官方发行版,根据薛定宇教授在其几本关于MA TLAB的著作中提到的,该函数返回的sol中结构体sol.x 和sol.y均为按行排列,与ode45等不同,不太规范(没办法,因为这个函数本来就不是Mathworks的官方作品),不过这一点已经不大可能得到改进了,因为L.F. Shampine 和S. Thompson 已经决定停止维护这个文件。
如果您想进一步了解该函数,可以访问它的主页。
MA TLAB帮助中关于该函数的介绍不很清楚,如果需要进一步了解这个函数,需要下载作者为其写的手册。
下面以MA TLAB 中所附的一个例程来说明这个函数与Simulink建模求解的不同。
在MA TLAB prompt中键入edit ddex1就会找看到函数作者所写的一个入门例子:function ddex1%DDEX1 Example 1 for DDE23.% This is a simple example of Wille' and Baker that illustrates the% straightforward formulation, computation, and plotting of the solution% of a system of delay differential equations (DDEs).%% The differential equations%% y'_1(t) = y_1(t-1)% y'_2(t) = y_1(t-1)+y_2(t-0.2)% y'_3(t) = y_2(t)%% are solved on [0, 5] with history y_1(t) = 1, y_2(t) = 1, y_3(t) = 1 for% t <= 0.%% The lags are specified as a vector [1, 0.2], the delay differential% equations are coded in the subfunction DDEX1DE, and the history is% evaluated by the function DDEX1HIST. Because the history is constant it% could be supplied as a vector:% sol = dde23(@ddex1de,[1, 0.2],ones(3,1),[0, 5]);%% See also DDE23, FUNCTION_HANDLE.% Jacek Kierzenka, Lawrence F. Shampine and Skip Thompson% Copyright 1984-2004 The MathW orks, Inc.% $Revision: 1.2.4.2 $ $Date: 2005/06/21 19:24:16 $sol = dde23(@ddex1de,[1, 0.2],@ddex1hist,[0, 5]);figure;plot(sol.x,sol.y)title('An example of Wille'' and Baker.');xlabel('time t');ylabel('solution y');% --------------------------------------------------------------------------function s = ddex1hist(t)% Constant history function for DDEX1.s = ones(3,1);% --------------------------------------------------------------------------function dydt = ddex1de(t,y,Z)% Differential equations function for DDEX1.ylag1 = Z(:,1);ylag2 = Z(:,2);dydt = [ ylag1(1)ylag1(1) + ylag2(2)y(2) ];这里先不管函数使用的具体语法,求解模型为:显然有两个延时常数1、0.2。
龙格-库塔法MATLAB

1. matlab 新建.m文件,编写龙格-库塔法求解函数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;%按照龙格库塔方法进行数值求解end2.另外再新建一个.,m文件,定义要求解的常微分方程函数function dx=fun1(t,x)dx =zeros(2,1);%初始化列向量dx(1) =0.08574*x(2)-1.8874*x(1)-20.17;dx(2) =1.8874*x(1)-0.08574*x(2)+115.16;3,再新建一个.m文件,利用龙格-库塔方法求解常微分方程[T1,F1]=runge_kutta1(@fun1,[46.30 1296 ],1,0,20); %求解步骤2定义的fun1常微分方程,@fun1是调用其函数指针,从0到20,间隔为1subplot(122)plot(T1,F1)%自编的龙格库塔函数效果title('自编的龙格库塔函数')grid运行步骤3文件即可得到结果,F1为估测值或者可以调用matlab自带函数ode45求解方法如下:[T,F] = ode45(@fun1,0:1:20,[17.64 37800 ]); subplot(121)plot(T,F)%Matlab自带的ode45函数效果title('ode45函数效果')grid。
matlab中三元二阶微分方程ode45编程

matlab中三元二阶微分方程ode45编程在Matlab中,可以使用ode45函数来求解三元二阶微分方程。
ode45是一个常用的数值解微分方程的函数,它使用了一种基于龙格-库塔方法的算法来求解微分方程。
下面是使用ode45函数求解三元二阶微分方程的步骤:1. 定义一个匿名函数,该函数输入参数为t和y,其中t表示时间,y 表示状态变量向量。
该匿名函数需要返回状态变量向量y的一阶导数dy/dt和二阶导数d2y/dt2。
例如,假设要求解以下三元二阶微分方程:d2x/dt2 + 2*dx/dt + 3*x = sin(t)则可以定义如下匿名函数:f = @(t,y) [y(2); -2*y(2)-3*y(1)+sin(t)];其中,y(1)表示状态变量x,y(2)表示状态变量dx/dt。
2. 定义初始条件。
即定义状态变量向量在初始时刻的值。
例如:y0 = [0;0];表示在初始时刻,状态变量x和dx/dt均为0。
3. 定义时间区间。
即定义需要求解微分方程的时间区间。
例如:tspan = [0,10];表示需要求解从时间0到时间10的微分方程。
4. 调用ode45函数进行求解。
例如:[t,y] = ode45(f,tspan,y0);其中,t表示求解得到的时间向量,y表示求解得到的状态变量向量矩阵。
矩阵y的每一行对应一个时间点的状态变量值。
5. 绘制结果。
可以使用plot函数绘制状态变量随时间的变化曲线。
例如:plot(t,y(:,1),'-o',t,y(:,2),'-*');表示绘制状态变量x和dx/dt随时间的变化曲线。
以上就是使用ode45函数求解三元二阶微分方程的基本步骤。
需要注意的是,当定义匿名函数时,需要按照正确的语法规则编写微分方程,并确保输入参数和输出结果符合要求。
此外,在定义初始条件和时间区间时也需要注意输入参数的格式和范围。
总之,使用ode45函数求解三元二阶微分方程是一种简单而有效的数值解法,可以帮助工程师和科学家快速得到微分方程的数值解,并用于实际应用中。
(完整word版)龙格库塔法求微分方程matlab

(完整word 版)龙格库塔法求微分方程matlab龙格—库塔方法求解微分方程初值问题(数学1201+41262022+陈晓云)初值问题:y x x -+=2dxdy ,10≤≤x 1)0(y = 四阶龙格-库塔公式:()y x K n n ,f 1=⎪⎪⎪⎪⎭⎫ ⎝⎛+=+K h y x K n h n 122f ,2 ⎪⎭⎫ ⎝⎛++=K y x fK h n h n 232,2 ()K h y h x f K n n 34,++=()K K K K y y h n 43211n 226++++=+ 程序:1)建立四阶龙格-库塔函数function [ x,y ] = nark4( dyfun,xspan,y0,h )% dyfun 为一阶微分方程的函数;y0为初始条件;xspan 表示x 的区间;h 为区间的步长; x=xspan(1):h:xspan(2);y(1)=y0;for n=1:length(x)-1k1=feval(dyfun,x(n),y(n));k2=feval(dyfun,x(n)+h/2,y(n)+h/2*k1);(完整word版)龙格库塔法求微分方程matlab 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*(k1+k2*2+2*k3+k4)/6;endx=x;y=y;2)执行程序(m文件)dyfun=inline('x^2+x-y');[x,y1]=nark4(dyfun,[0,1],1,0.1);x=0:0.1:1;Format longy2=x.^2-x+1R4=y2-y1[x',y1',y2',R4']y2=dsolve('Dy=x^2+x-y','y(0)=1','x')plot(x,y1,'b*-')hold ony3=inline('x^2-x+1')fplot(y3,[0,1],'ro-')legend('R-K4','解析解')3)执行结果ans =X RK4近似值解析值0 1.000000000000000 1.000000000000000 0.100000000000000 0.910000208333333 0.910000000000000 0.200000000000000 0.840000396841146 0.840000000000000 0.300000000000000 0.790000567410084 0.790000000000000 0.400000000000000 0.760000721747255 0.760000000000000 0.500000000000000 0.750000861397315 0.750000000000000 0.600000000000000 0.760000987757926 0.760000000000000 0.700000000000000 0.790001102093746 0.790000000000000 0.800000000000000 0.840001205549083 0.8400000000000000.900000000000000 0.910001299159352 0.9100000000000001.000000000000000 1.000001383861433 1.000000000000000误差-0.000000208333333-0.000000396841146-0.000000567410084-0.000000721747255-0.000000861397315-0.000000987757926-0.000001102093746-0.000001205549083-0.000001299159352-0.000001383861433 y2 =x^2 - x + 1结果分析:初值问题的解析解为Y=x^2 - x + 1由图看出龙格库塔方法误差很小,具有很高的精度。
matlab中常用来解常微分方程的原理

matlab中常用来解常微分方程的原理Matlab是一种常用的科学计算软件,它在解常微分方程方面具有很强的功能。
在实际应用中,常微分方程是描述自然现象和工程问题的重要数学模型。
通过使用Matlab,我们可以方便地求解各种类型的常微分方程,从而得到问题的解析解或数值解。
本文将介绍Matlab中常用来解常微分方程的原理和方法。
常微分方程是描述一个未知函数及其导数之间关系的方程。
它可以分为初值问题和边值问题两种类型。
初值问题是给定方程的初始条件,要求求解出函数的解析解或数值解;而边值问题是给定方程在两个或多个点上的边界条件,要求求解出满足这些条件的函数解。
Matlab提供了多种方法来解这些问题,如常微分方程求解器、符号计算工具箱和数值计算工具箱等。
我们来介绍常微分方程求解器。
Matlab中的常微分方程求解器可以直接求解一阶或高阶的常微分方程。
常见的求解器有ode45、ode23、ode15s等。
其中,ode45是一种常用的求解器,它采用的是龙格-库塔方法,可以求解刚性和非刚性的常微分方程。
使用这些求解器,我们只需要定义好方程的形式和初始条件,就可以得到方程的数值解。
我们介绍符号计算工具箱。
符号计算工具箱可以对常微分方程进行符号计算,得到方程的解析解。
通过符号计算,我们可以得到方程的精确解,而不仅仅是数值解。
在Matlab中,符号计算工具箱提供了dsolve函数来求解常微分方程。
我们只需要输入方程的形式,就可以得到方程的解析解。
符号计算工具箱对于一些简单的常微分方程求解非常方便,但对于复杂的方程,可能需要更复杂的方法来求解。
我们介绍数值计算工具箱。
数值计算工具箱提供了各种数值方法来求解常微分方程。
常见的数值方法有欧拉方法、龙格-库塔方法和有限差分法等。
这些方法可以将常微分方程转化为一系列迭代计算,从而得到方程的数值解。
在Matlab中,数值计算工具箱提供了多个函数来实现这些数值方法,如euler、ode23s和fdcoeffF等。
【免费下载】数值分析Matlab作业龙格库塔欧拉方法解二阶微分方程

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('欧拉');
稳定性很好,RK4 法是四阶方法,每步的误差是 h5 阶,而总积累误差为 h4 阶。所以比欧 拉稳定。
运行第三个程序:在一幅图中显示欧拉法和 RK4 法,随着截断误差的积累,欧拉法产生了 较大的误差 h=0.01
对全部高中资料试卷电气设备,在安装过程中以及安装结束后进行高中资料试卷调整试验;通电检查所有设备高中资料电试力卷保相护互装作置用调与试相技互术关,系电,力根通保据过护生管高产线中工敷资艺设料高技试中术卷资,配料不置试仅技卷可术要以是求解指,决机对吊组电顶在气层进设配行备置继进不电行规保空范护载高与中带资负料荷试下卷高问总中题体资,配料而置试且时卷可,调保需控障要试各在验类最;管大对路限设习度备题内进到来行位确调。保整在机使管组其路高在敷中正设资常过料工程试况中卷下,安与要全过加,度强并工看且作护尽下关可都于能可管地以路缩正高小常中故工资障作料高;试中对卷资于连料继接试电管卷保口破护处坏进理范行高围整中,核资或对料者定试对值卷某,弯些审扁异核度常与固高校定中对盒资图位料纸置试,.卷保编工护写况层复进防杂行腐设自跨备动接与处地装理线置,弯高尤曲中其半资要径料避标试免高卷错等调误,试高要方中求案资技,料术编试交写5、卷底重电保。要气护管设设装线备备置敷4高、调动设中电试作技资气高,术料课中并3中试、件资且包卷管中料拒含试路调试绝线验敷试卷动槽方设技作、案技术,管以术来架及避等系免多统不项启必方动要式方高,案中为;资解对料决整试高套卷中启突语动然文过停电程机气中。课高因件中此中资,管料电壁试力薄卷高、电中接气资口设料不备试严进卷等行保问调护题试装,工置合作调理并试利且技用进术管行,线过要敷关求设运电技行力术高保。中护线资装缆料置敷试做设卷到原技准则术确:指灵在导活分。。线对对盒于于处调差,试动当过保不程护同中装电高置压中高回资中路料资交试料叉卷试时技卷,术调应问试采题技用,术金作是属为指隔调发板试电进人机行员一隔,变开需压处要器理在组;事在同前发一掌生线握内槽图部内 纸故,资障强料时电、,回设需路备要须制进同造行时厂外切家部断出电习具源题高高电中中源资资,料料线试试缆卷卷敷试切设验除完报从毕告而,与采要相用进关高行技中检术资查资料和料试检,卷测并主处且要理了保。解护现装场置设。备高中资料试卷布置情况与有关高中资料试卷电气系统接线等情况,然后根据规范与规程规定,制定设备调试高中资料试卷方案。
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*的近似值,显然能构造出具有很高精度的高阶计算公式。
matlab龙格库塔法解微分方程组

matlab龙格库塔法解微分方程组
一、引言
龙格库塔法是数值计算中常用的一种求解微分方程的方法,其具有较高的精度和稳定性。
在MATLAB中,可以使用ode45函数来实现龙格库塔法求解微分方程组。
二、龙格库塔法简介
龙格库塔法是一种常用的数值积分方法,也可用于求解微分方程。
该方法将微分方程转化为一个初值问题,并采用逐步逼近的方式计算出数值解。
三、使用ode45函数求解微分方程组
在MATLAB中,可以使用ode45函数来求解微分方程组。
该函数使用了龙格库塔法进行数值计算,并提供了较高的精度和稳定性。
四、MATLAB代码实现
以下是一个使用ode45函数求解微分方程组的示例代码:
function dydt = myfun(t,y)
dydt = zeros(2,1);
dydt(1) = y(2);
dydt(2) = -sin(y(1));
end
[t,y] = ode45(@myfun,[0 10],[0;1]);
plot(t,y(:,1),'-o',t,y(:,2),'-o')
xlabel('Time')
ylabel('Solution')
legend('y_1','y_2')
五、总结
龙格库塔法是一种常用的数值计算方法,可以用于求解微分方程。
在MATLAB中,可以使用ode45函数来实现龙格库塔法求解微分方程组。
通过以上示例代码,我们可以看到MATLAB提供了较为简单的方式来实现龙格库塔法求解微分方程组,并且具有较高的精度和稳定性。
龙格-库塔法,求解常微分方程

隆格库塔法求解常微分方程摘要科学技术中常常需要求解常微分方程的定解问题,这里问题最简单的形式,是本章将着重考察的一阶方程的初值问题.虽然求解常微分方程有各种各样的解析方法,但解析方法只能用来求解一些特殊类型的方程,实际问题中归结出来的微分方程主要靠数值解法求解.本文着重讨论了隆格库塔法求解一阶常微分方程的初值问题,采用了精度较高的经典的四阶隆格库塔法,然后通过对实例运用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 应用使用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');。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
龙格-库塔方法是一种经典方法,具有很高的精度,它间接的利用了泰勒级数展开,避免了高阶偏导数的计算。
此处以最为经典的四级四阶龙格-库塔方法为例,计算格式如下
()()()11234121324
3226,,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 +⎧=++++⎪⎪⎪=⎪⎪⎛⎫=++⎨ ⎪⎝⎭⎪⎪⎛⎫=+⎪ ⎪⎝⎭⎪⎪=++⎩
1龙格-库塔法解一阶ODE 对于形如()()0, dy f x y a x b dx y a y ⎧=<≤⎪⎨⎪=⎩
的一阶ODE 初值问题,可以直接套用公式,如今可以借助计算机方便的进行计算,下面给出一个实例
()2 0101dy x y x dx y y ⎧=-<≤⎪⎨⎪=⎩
取步长h=0.1
,此处由数学知识可得该方程的精确解为y =。
在这里利用MATLAB 编程,计算数值解并与精确解相比,代码如下:
(1)写出微分方程,便于调用和修改
function val = odefun( x,y )
val = y-2*x/y;
end
(2)编写runge-kutta 方法的函数代码
function y = runge_kutta( h,x0,y0 )
k1 = odefun(x0,y0);
k2 = odefun(x0+h/2,y0+h/2*k1);
k3 = odefun(x0+h/2,y0+h/2*k2);
k4 = odefun(x0+h,y0+h*k3);
y = y0+h*(k1+2*k2+2*k3+k4)/6;
end
(3)编写主函数解微分方程,并观察数值解与精确解的差异clear all
h = 0.1;
x0 = 0;
y0 = 1;
x = 0.1:h:1;
y(1) = runge_kutta(h,x0,y0);
for k=1:length(x)
x(k) = x0+k*h;
y(k+1) = runge_kutta(h,x(k),y(k));
end
z = sqrt(1+2*x);
plot(x,y,’*’);
hold on
plot(x,z,'r');
结果如下图,数值解与解析解高度一致
2龙格-库塔法解高阶ODE
对于高阶ODE来说,通用的方法是将高阶方程通过引入新的变量降阶为一阶方程组,此处仍以一个实例进行说明。
++=
5002007502000
y y y
这是一个二阶ODE,描述的是一个物体的有阻尼振动情况。
初始条件为000 0y y == ,将方程降阶,引入一个向量型变量Y
y Y y ⎡⎤=⎢⎥⎣⎦故有()2000200750500y y dY Y y y y dt ⎡⎤⎡⎤⎢⎥===--⎢⎥⎢⎥⎣⎦⎢⎥⎣⎦
记()()1 2y Y y Y ==则()()()()2200020027501500Y Y Y Y ⎡⎤⎢⎥=--⎢⎥⎢⎥⎣⎦
至此,二阶方程降阶为一阶
方程组。
值得注意的是此时再用龙格-库塔法进行求解时,代入的将是一个Y 向量。
同样利用MATLAB 进行计算,步长h=0.05,时间周期为[0,20].
(1) 编写ODE 函数
function Y = odefun1( ~,Y0 )
% 此处Y0为一个列向量,因为时间t 未显含在一阶方程组中 % 所以ode 函数的第一个参数为空,要根据具体情况而定。
Y = [Y0(2);
(2000-200*Y0(2)-750*Y0(1))/500;];
end
(2) 编写runge-kutta 函数
function Y = rkfa( h,t0,Y0 )
k1 = odefun1(t0,Y0);
k2 = odefun1(t0+h/2,Y0+h/2*k1);
k3 = odefun1(t0+h/2,Y0+h/2*k2);
k4 = odefun1(t0+h,Y0+h*k3);
Y = Y0+h*(k1+2*k2+2*k3+k4)/6;
end
(3)编写主函数
clear all
h = 0.05;
t = 0.05:h:20;
t0 = 0;
Y0 = [0;
0];%初值
Y = cell(1,length(t));
Y{1} = rkfa( h,t0,Y0 );
z = zeros(2,length(t));
for k=1:length(t)
Y{k+1} = rkfa( h,t0,Y{k}); z(1,k) = Y{k}(1);
z(2,k) = Y{k}(2);
end
plot(t,z(1,:),'r');%位移y的图像hold on
plot(t,z(2,:));%速度y’的图像
求解结果如下图。