matlab编的4阶龙格库塔法解微分方程的程序
MATLAB改进欧拉法与四阶龙格-库塔求解一阶常微分方程
姓名:樊元君学号:02 日期:
一、实验目的
掌握MATLAB语言、C/C++语言编写计算程序的方法、掌握改进欧拉法与四阶龙格-库塔求解一阶常微分方程的初值问题。掌握使用MATLAB程序求解常微分方程问题的方法。
:
二、实验内容
1、分别写出改进欧拉法与四阶龙格-库塔求解的算法,编写程序上机调试出结果,要求所编程序适用于任何一阶常微分方程的数值解问题,即能解决这一类问题,而不是某一个问题。
实验中以下列数据验证程序的正确性。
求,步长h=。
*
2、实验注意事项
的精确解为,通过调整步长,观察结果的精度的变化
^
)
三、程序流程图:
●改进欧拉格式流程图:
~
|
●四阶龙格库塔流程图:
]
四、源程序:
●改进后欧拉格式程序源代码:
function [] = GJOL(h,x0,y0,X,Y)
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);
使用Matlab进行微分方程求解的方法
使用Matlab进行微分方程求解的方法引言
微分方程是数学中非常重要的一部分,广泛应用于物理、经济、工程等领域。
对于大部分微分方程的解析解往往难以求得,而数值解法则成为了一种常用的解决手段。Matlab作为一种强大的科学计算软件,也提供了丰富的工具和函数用于求
解微分方程,本文将介绍一些常见的使用Matlab进行微分方程求解的方法。
一、数值求解方法
1. 欧拉方法
欧拉方法是最简单的一种数值求解微分方程的方法,它将微分方程的微分项用
差分的方式进行近似。具体的公式为:
y(n+1) = y(n) + hf(x(n), y(n))
其中,y(n)表示近似解在第n个点的值,h为步长,f(x, y)为微分方程的右端项。在Matlab中使用欧拉方法进行求解可以使用ode113函数,通过设定不同的步长,
可以得到不同精度的数值解。
2. 中点法
中点法是较为精确的一种数值求解微分方程的方法,它的计算公式为:
k1 = hf(x(n), y(n))
k2 = hf(x(n) + h/2, y(n) + k1/2)
y(n+1) = y(n) + k2
中点法通过计算两个斜率的平均值来得到下一个点的值,相较于欧拉方法,中
点法能提供更精确的数值解。
3. 4阶龙格库塔法
龙格库塔法是一类高阶数值求解微分方程的方法,其中4阶龙格库塔法是最常
用的一种。它的计算公式为:
k1 = hf(x(n), y(n))
k2 = hf(x(n) + h/2, y(n) + k1/2)
k3 = hf(x(n) + h/2, y(n) + k2/2)
k4 = hf(x(n) + h, y(n) + k3)
MATLAB实验四_求微分方程的解
ode45
非刚性
单步法;4,5 阶 R-K 方法;大部分场合的首选方法 累计截断误差为 (△x)3
ode23
ode113
非刚性
非刚性
单步法;2,3 阶 R-K 方法;使用于精度较低的情形 累计截断误差为 (△x)3
多步法;Adams算法;高低精 计算时间比 ode45 短 度均可到 10-3~10-6
n
xn1 xn b k 0, 1, 2, , n 1
步长:h xk 1 xk (b a) / n, 差商代替微商
y( xk 1 ) y( xk ) dy y( xk 1 ) y( xk ) h y '( xk ) dx x h k k = 0, 1, 2, ..., n-1 y0 y( x0 ) 得方程组: yk 1 yk h f ( xk , yk ) x x h y 是 y (x ) 的近似 k k 1
dsolve('Dy=2*x','x'); dsolve('Dy=2*x');
t
% dy/dx = 2x % dy/dt = 2x
若找不到解析解,则返回其积分形式。
dsolve 举例
例 2:求微分方程 xy ' y e x 0 在初值条件 y(1) 2e
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/h
x=(i-1)*h;
K=f(x,y);
y=y+h*K;
fprintf('欧拉法:x=%f, y=%f\n',x+h,y);
end
fprintf('\n');
%改进的欧拉法
y=1;
fprintf('改进的欧拉法:x=%f, y=%f\n',0,1);
for i=1:1/h
x=(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);
end
fprintf('\n');
%龙格库塔⽅法
y=1;
fprintf('龙格库塔法:x=%f, y=%f\n',0,1);
for i=1:1/h
x=(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);
四阶龙格库塔法(Runge-Kutta)求解微分方程
四阶龙格库塔法(Runge-Kutta )求解微分方程
张晓颖
(天津大学 材料学院 学号:1012208027)
1 引言
计算传热学中通常需要求解常微分方程。这类问题的简单形式如下:
{
),(')(0
0y x f y y x y == (1)
虽然求解常微分方程有各种各样的解析方法,但解析方法只能用来求解一些特殊类型的方程,实际问题中的多数微分方程需要采用数值解法求解。初值问题(1)的数值解法有个基本特点,它们采取“步进式”,即求解过程顺着节点排序一步一步向前推进。这类算法是要给出用已知信息y n 、 y n −i ……计算y n +1的递推公式即可。
2 龙格库塔法(Runge-Kutta )介绍
假设对于初值问题(1)有解 y = y (x ) ,用 Taylor 展开有:
......)(!
3)(!2)()()(3
21+'''+''+'+=+n n n n n x y h x y h x y h x y x y (2) 龙格库塔法(Runge-Kutta )实质上是间接的使用 Taylor 级数法的一种方法。对于差商
h
x y x y n n )
()(1-+,根据微分中值定理,存在 0 < θ < 1 ,使得:
)()
()(1h x y h
x y x y n n θ+'=-+ (3)
于是对于 y = y (x ) ,可得到:
))(,()()(1h x y h x hf x y x y n n n n θθ+++=+ (4)
设))(,(*
h x y h x f K n n θθ++=,为区间 [x n , x n +1 ] 上的平均斜率。四阶龙
龙格库塔法求微方程matlab
龙格—库塔方法求解微分方程初值问题
(数学1201+41262022+陈晓云)
初值问题:
y x x -+=2dx dy ,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 f K h n h n 232,2
()K h y h x f K n n 34,++=
()K K K K y y h n 4
3211n 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)-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*(k1+k2*2+2*k3+k4)/6;
end
x=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 long
matlab4阶龙格库塔法程序
====Word行业资料分享--可编辑版本--双击可删====
function varargout=saxplaxliu(varargin)
clc,clear
x0=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:n
fprintf('%2d %4.4f %4.4f\n',i,x(i),y(i));
end
function 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-1
K1=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);
end
y=y1;
结果:
源-于-网-络-收-集
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函数,可以快速求解常微分方程的数值解,并通过绘图来展示结果。希望本篇文章对读者
matlab4阶龙格库塔法程序
function varargout=saxplaxliu(varargin)
clc,clear
x0=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:n
fprintf('%2d %4.4f %4.4f\n',i,x(i),y(i)); end
function 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-1
K1=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);
end
y=y1;
结果:
i x(i) y(i)
1 0.0000 1.0000
2 0.1000 0.9901
3 0.2000 0.9615
4 0.3000 0.9174
5 0.4000 0.8621
6 0.5000 0.8000
7 0.6000 0.7353
8 0.7000 0.6711
9 0.8000 0.6098
10 0.9000 0.5525
11 1.0000 0.5000
12 1.1000 0.4525
13 1.2000 0.4098
四阶龙格库塔法求解微分方程组
四阶龙格库塔法求解微分方程组
四阶龙格库塔法是一种数值解微分方程组的方法。它的基本思想
是将微分方程组中的每个方程都离散化,并使用预测-校正的手段来逐
步计算方程组的数值解。
具体地,四阶龙格库塔法可以分为以下步骤:
1. 对微分方程组进行离散化,即将其转化为一组差分方程。这
一步骤的详细方法因具体的微分方程组而异,需要根据情况进行逐步
求解。
2. 使用四个预测基于中间值的校正步骤,来计算方程组的数值解。具体地,假设当前时刻为t,并已知方程组在时刻t的解为y(t),则可以按照以下步骤逐步计算y(t+Δt):
(1)预测y(t+Δt)的初值为y1 = y(t) + (Δt/2)k1
(2)计算y1对应的导数k2 = f(t+Δt/2, y1)
(3)校正y(t+Δt)的初值为y2 = y(t) + (Δt/2)k2
(4)计算y2对应的导数k3 = f(t+Δt/2, y2)
(5)再次校正y(t+Δt)的初值为y3 = y(t) + Δt k3
(6)计算y3对应的导数k4 = f(t+Δt, y3)
(7)利用k1, k2, k3, k4来计算y(t+Δt),即y(t+Δt) = y(t) + (Δt/6)(k1+2k2+2k3+k4)
3. 重复第二步,直到计算出方程组在终止时刻t1的解y(t1)。
重要的是控制步长,以避免数值计算精度不足或者计算速度过慢。
综上所述,四阶龙格库塔法是一种比较常用的解微分方程组的数
值算法。但由于龙格库塔法的计算步骤较多,且需要进行一定的离散
化处理,因此具有一定的计算复杂度和技巧难度。
四阶龙格库塔法求解微分方程组
四阶龙格库塔法求解微分方程组
微分方程是数学中的一个重要分支,它描述了自然界中许多现象的发展规律。在实际应用中,我们经常需要求解微分方程组,以得到系统的演化规律和性质。本文将介绍一种常用的求解微分方程组的数值方法——四阶龙格库塔法。
一、微分方程组的数值解法
微分方程组是形如下式的方程集合:
begin{cases} frac{dy_1}{dx}=f_1(x,y_1,y_2,cdots,y_n) frac{dy_2}{dx}=f_2(x,y_1,y_2,cdots,y_n) cdots
frac{dy_n}{dx}=f_n(x,y_1,y_2,cdots,y_n) end{cases} 其中,$y_1,y_2,cdots,y_n$是关于自变量$x$的未知函数,
$f_1,f_2,cdots,f_n$是已知的函数。求解微分方程组的数值方法主要有以下两种:
1.欧拉法
欧拉法是最简单的数值方法之一,它的基本思想是将微分方程组中的导数用差分代替,从而得到一个递推公式。具体而言,欧拉法的递推公式为:
$$y_{i+1}=y_i+hf(x_i,y_i)$$
其中,$h$是步长,$x_i$和$y_i$分别是第$i$个点的自变量和因变量的值。欧拉法的精度较低,容易出现数值误差,但是它的计算速度快,适用于一些简单的微分方程组。
2.龙格库塔法
龙格库塔法是一种常用的高精度数值方法,它的基本思想是将微分方程组中的导数用一系列加权的差分代替,从而得到一个递推公式。其中,四阶龙格库塔法是最为常用的一种。具体而言,四阶龙格库塔法的递推公式为:
龙格-库塔法(Runge-Kutta)matlab代码及含义
龙格-库塔法(Runge-Kutta)
数值分析中,龙格-库塔法(Runge-Kutta)是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。这些技术由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。
经典四阶龙格库塔法
RK4””或者就是龙格库塔法的家族中的一个成员如此常用,以至于经常被称为“RK4
“龙格库塔法”。
令初值问题表述如下。
则,对于该问题的RK4由如下方程给出:
其中
这样,下一个值(yn+1)由现在的值(yn)加上时间间隔(h)和一个估算的斜率的乘积决定。该斜率是以下斜率的加权平均:
k1是时间段开始时的斜率;
k2是时间段中点的斜率,通过欧拉法采用斜率k1来决定y在点tn+h/2的值;k3也是中点的斜率,但是这次采用斜率k2决定y值;
k4是时间段终点的斜率,其y值用k3决定。
当四个斜率取平均时,中点的斜率有更大的权值:
RK4法是四阶方法,也就是说每步的误差是h5阶,而总积累误差为h4阶。
注意上述公式对于标量或者向量函数(y可以是向量)都适用。
显式龙格库塔法
显示龙格-库塔法是上述RK4法的一个推广。它由下式给出
其中
(注意:上述方程在不同著述中由不同但却等价的定义)。
要给定一个特定的方法,必须提供整数s(阶段数),以及系数aij(对于1≤j
c3a31
如果要求方法有精度p则还有相应的条件,也就是要求舍入误差为O(hp+1)时的条件。这些可以从舍入误差本身的定义中导出。例如,一个2阶精度的2段方法要求b1+b2=1,b2c2=1/2,以及b2a21=1/2。
经典四阶龙格库塔法解一阶微分方程组
数值计算课程设计
1. 经典四阶龙格库塔法解一阶微分方程组
1.1 运用四阶龙格库塔法解一阶微分方程组算法分析
x k 1 x k h
(f 1 2 f 2 2 f 3 f 4)
6
h
y k 1 y k h 6
(g 1 2g 2 2g 3 g 4 )
t k 1 t k h
经过循环计算由 t 0
,x 0
, y 0
推得 t 1
,x 1
,y
1 t 2
,x 2,y 2
⋯⋯
每个龙格 -库塔方法都是由一个合适的泰勒方法推导而来,使得其最终全局 误差为 O h N
, 一种折中方法是每次进行若干次函数求值,从而省去高阶导
数计 算。4阶龙格-库塔方法 (RK4)是最常用的,它适用于一般的应用,因为它非常精 准,稳定,且易于编程。
f 1 f (t k , x k , y k ) ,
h
f 2 f (t k h 2
,x k
hh
2
f 1,y k 2
g 1)
h
f 3 f (t k h 2
,x k
h
f 2,y k 2
g 2)
f 4 f (
t k h ,
x k hf 3,
y k hg 3)
g 1 g(t k ,x k , y k ) h g 2
g(t k 2
,x k
hh
2
f 1,y k 2
g 1)
h
g 3 g(t k 2
, x k
h
f 2
h
2,y k 2
g 2) g 4 g(t k h,x k hf 3, y k hg 3)
1-1)
1-2)
1-3)
1-4)
1-5)
1-6)
1-7)
1-8)
1-9)
1-10 )
经典四阶龙格库塔法解一阶微分方程组
1.2 经典四阶龙格库塔法解一阶微分方程流程图
1.3 经典四阶龙格库塔法解一阶微分方程程序代码:
四阶龙格库塔的MATLAB实现
一阶微分方程的四阶龙格库塔的MATLAB实现
作者:头铁的小甘
当计算机求解常系数微分方程时,一般存在龙格库塔法和有限元分析这两种方法。龙格库塔法基于微分中值定理,而有限元差分法基于泰勒展开。本篇文章实现MATLAB实现龙格库塔解法。简单公式如下,(详细内容自己百度或者其他来源)
ⅆy
ⅆx
=f(x,y) y(0)=y0
y j+1=y j+1
6
(k1+2k2+2k3+k4)ℎk1=f(x i,y i)
k2=f(x i+1
ℎ,y i+
1
k1ℎ)
k3=f(x i+1
2
ℎ,y i+
1
2
k2ℎ)
k4=f(x i+ℎ,y i+k3ℎ)
其中h为步长
物理问题背景
现在我们假设解决一个简单的RC串联电路,输出选取电容电容两端的电压,
输入电压选取方波,电路图如下所示:
微分方程为
ⅆU c ⅆt =
U i−U c
RC
其中U c为电容两端的电压, U I为输入电压。MATLAB实现满足下面流程图
MATLAB代码
图1
图2
图3
图1为主程序,图2为龙格库塔法,图3为微分方程。
结果显示
下一篇:二阶常微分方程的龙格库塔法求解及其MATLAB实现
matlab编的4阶龙格库塔法解微分方程的程序
matlab编的4阶龙格库塔法解微分方程的程序
2010-03—10 20:16
function varargout=saxplaxliu(varargin)
clc,clear
x0=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:n
fprintf('%2d %4.4f %4.4f\n’,i,x(i),y(i));end
function 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-1
K1=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);
end
y=y1;
结果:
i x(i) y(i)
1 0.0000 1。0000
2 0。1000 0.9901
3 0。2000 0.9615
4 0.3000 0。9174
5 0.4000 0。8621
6 0.5000 0。8000
7 0。6000 0。7353
8 0。7000 0.6711
9 0。8000 0。6098
10 0。9000 0。5525
11 1.0000 0.5000
matlab用四阶龙格库塔函数求解微分方程组
一、介绍
Matlab作为一种强大的科学计算软件,提供了众多函数和工具来解决微分方程组。其中,四阶龙格库塔函数是一种常用的数值方法,用于求解常微分方程组。本文将介绍如何使用Matlab中的四阶龙格库塔函数来求解微分方程组,并对该方法的原理和实现进行详细说明。
二、四阶龙格库塔方法
四阶龙格库塔方法是一种常用的数值方法,用于求解常微分方程组。它是一种显式的Runge-Kutta方法,通过逐步逼近微分方程的解,在每一步使用多个中间值来计算下一步的解。该方法通过四个中间值来计算下一步的状态,并且具有较高的精度和稳定性。
三、在Matlab中使用四阶龙格库塔方法求解微分方程组
在Matlab中,可以使用ode45函数来调用四阶龙格库塔方法来解决微分方程组的问题。ode45函数是Matlab提供的用于求解常微分方程组的函数,可以通过指定微分方程组以及初值条件来调用四阶龙格库塔方法来进行求解。
1. 定义微分方程组
我们需要定义要求解的微分方程组。可以使用Matlab中的匿名函数来定义微分方程组,例如:
```matlab
f = (t, y) [y(2); -sin(y(1))];
```
其中,f是一个匿名函数,用于表示微分方程组。在这个例子中,微分方程组是y' = y2, y2' = -sin(y1)。
2. 指定初值条件和求解区间
接下来,我们需要指定微分方程组的初值条件和求解区间。初值条件可以通过指定一个初始时刻的状态向量来完成,例如:
```matlab
tspan = [0, 10];
y0 = [0, 1];
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
matlab编的4阶龙格库塔法解微分方程的程序
2010-03-10 20:16
function varargout=saxplaxliu(varargin)
clc,clear
x0=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:n
fprintf('%2d %4.4f %4.4f\n',i,x(i),y(i)); end
function 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-1
K1=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);
end
y=y1;
结果:
i x(i) y(i)
1 0.0000 1.0000
2 0.1000 0.9901
3 0.2000 0.9615
4 0.3000 0.9174
5 0.4000 0.8621
6 0.5000 0.8000
7 0.6000 0.7353
8 0.7000 0.6711
9 0.8000 0.6098
10 0.9000 0.5525
11 1.0000 0.5000
12 1.1000 0.4525
13 1.2000 0.4098
龙格库塔法
一、基本原理:
龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高,采取措施对误差进行抑制,所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。对于一阶精度的欧拉公式有:
yi+1=yi+h*K1
K1=f(xi,yi)
当用点xi处的斜率近似值K1与右端点xi+1处的斜率K2的算术平均值作为平均斜率K*的近似值,那么就会得到二阶精度的改进欧拉公式:
yi+1=yi+h*( K1+ K2)/2
K1=f(xi,yi)
K2=f(xi+h,yi+h*K1)
依次类推,如果在区间[xi,xi+1]内多预估几个点上的斜率值K1、K2、……Km,并用他们的加权平均数作为平均斜率K*的近似值,显然能构造出具有很高精度的高阶计算公式。经数学推导、求解,可以得出四阶龙格-库塔公式,也就是在工程中应用广泛的经典龙格-库塔算法:
yi+1=yi+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)
通常所说的龙格-库塔法是指四阶而言的,我们可以仿二阶、三阶的情形推导出常用的标准四阶龙格-库塔法公式
(1)
计算公式(1)的局部截断误差是。
龙格-库塔法具有精度高,收敛,稳定(在一定条件下),计算过程中可以改变步长,不需要计算高阶导数等优点,但仍需计算在一些点上的值,如四阶龙格-库塔法每计算一步需要计算四次的值,这给实际计算带来一定的复杂性,因此,多用来计算“表头”。
二、小程序
#include
#include
#define f(x,y) (-1*(x)*(y)*(y))
void main(void)
{
double a,b,x0,y0,k1,k2,k3,k4,h;
int n,i;
printf("input a,b,x0,y0,n:");
scanf("%lf%lf%lf%lf%d",&a,&b,&x0,&y0,&n);
printf("x0\ty0\tk1\tk2\tk3\tk4\n");
for(h=(b-a)/n,i=0;i!=n;i++)
{
k1=f(x0,y0);
k2=f(x0+h/2,y0+k1*h/2);
k3=f(x0+h/2,y0+k2*h/2);
k4=f(x0+h,y0+h*k3);
printf("%lf\t%lf\t",x0,y0);
printf("%lf\t%lf\t",k1,k2);
printf("%lf\t%lf\n",k3,k4);
y0+=h*(k1+2*k2+2*k3+k4)/6;
x0+=h;
}
printf("xn=%lf\tyn=%lf\n",x0,y0);
}
运行结果:
input a,b,x0,y0,n:0 5 0 2 20
x0 y0 k1 k2 k3 k4
0.000000 2.000000 -0.000000 -0.500000 -0.469238 -0.886131
0.250000 1.882308 -0.885771 -1.176945 -1.129082 -1.280060
0.500000 1.599896 -1.279834 -1.295851 -1.292250 -1.222728
0.750000 1.279948 -1.228700 -1.110102 -1.139515 -0.990162
1.000000 1.000027 -1.000054 -0.861368 -0.895837 -0.752852
1.250000 0.780556 -0.761584 -0.645858 -0.673410 -0.562189
1.500000 0.615459 -0.568185 -0.481668 -0.500993 -0.420537
1.750000 0.492374 -0.424257 -0.361915 -0.374868 -0.317855
2.000000 0.400054 -0.320087 -0.275466 -0.284067 -0.243598