矩阵龙格库塔

合集下载

龙格库塔法解微分方程组

龙格库塔法解微分方程组

龙格库塔法解微分方程组引言微分方程组是数学中经常遇到的问题,在物理、工程和自然科学中都有广泛的应用。

为了求解微分方程组,我们需要利用数值方法来逼近解析解。

本文将介绍一种常用的数值方法——龙格库塔法(Runge-Kutta method),并探讨如何利用该方法来解微分方程组。

龙格库塔法概述龙格库塔法是一种迭代数值方法,用于求解常微分方程的初值问题。

它的主要思想是将微分方程的解进行离散化,将其转化为一系列的逼近值。

龙格库塔法的基本步骤如下:1.确定步长h和迭代次数n。

2.初始化初始条件,并假设第一个逼近值为y(xi)。

3.依次计算每个逼近值,直到得到y(xi+n*h)为止。

在每个迭代步骤中,龙格库塔法根据前一步的逼近值来计算下一步的逼近值。

该方法具有高阶精度和较好的稳定性,在实际应用中广泛使用。

单一微分方程的龙格库塔法首先,我们来看如何利用龙格库塔法来解一阶常微分方程。

以方程dy/dx = f(x, y)为例,其中f(x, y)为给定的函数。

步骤一:确定步长和迭代次数选择合适的步长h和迭代次数n来进行数值计算。

步长h决定了离散化的精度,而迭代次数n决定了逼近解的数目。

步骤二:初始化条件并计算逼近值设初始条件为y(x0) = y0,其中x0为起始点,y0为起始点处的函数值。

我们先通过欧拉法计算出y(x0 + h)的逼近值,然后再通过该逼近值来计算下一个逼近值。

逼近值的计算公式如下:k1 = h * f(x0, y0)k2 = h * f(x0 + h/2, y0 + k1/2)k3 = h * f(x0 + h/2, y0 + k2/2)k4 = h * f(x0 + h, y0 + k3)y(x0 + h) = y0 + 1/6 * (k1 + 2k2 + 2k3 + k4)步骤三:重复步骤二直到得到y(xi+n*h)依次利用上一步计算出的逼近值来计算下一个逼近值,直到得到y(xi+n*h)为止。

微分方程组的龙格库塔法对于一阶微分方程组的初值问题,我们可以将其转化为向量形式。

龙格库塔法解轨道参数

龙格库塔法解轨道参数

龙格库塔法解轨道参数引言龙格库塔法(Runge-Kutta method)是一种常用的数值计算方法,用于求解常微分方程的数值解。

在天体力学中,我们经常需要通过数值方法来计算天体的轨道参数,如轨道椭圆的长短轴、离心率、倾角等。

本文将介绍龙格库塔法在解轨道参数中的应用,并详细探讨该方法的原理和实现过程。

基本原理龙格库塔法是一种迭代求解的方法,在每个时间步长内利用当前的状态来估计下一个状态。

具体而言,龙格库塔法将微分方程的求解问题转化为一个迭代的求解问题,通过逐步迭代来逼近精确解。

在解轨道参数的问题中,我们通常需要根据已知的初始条件以及天体的质量和力学模型来求解天体的轨道参数。

常用的力学模型有开普勒模型和牛顿模型。

龙格库塔法可以根据力学模型的不同进行相应的求解。

开普勒模型下的轨道参数求解步骤一:确定初始条件在使用龙格库塔法求解轨道参数之前,我们需要确定一些初始条件。

这些初始条件包括天体的质量、位置和速度。

步骤二:选择时间步长在求解过程中,我们需要选择一个合适的时间步长。

时间步长越小,计算的精度会越高,但计算的时间会增加。

步骤三:迭代求解利用龙格库塔法进行迭代求解的具体步骤如下:1.根据当前时刻的位置和速度,计算天体在该时刻的加速度。

2.根据当前时刻的位置、速度和加速度,计算下一个时刻的位置和速度。

3.更新当前时刻的位置和速度为新的位置和速度。

4.重复上述步骤,直到达到指定的终止条件。

步骤四:计算轨道参数通过迭代求解,我们可以得到天体在不同时刻的位置和速度。

根据这些位置和速度,我们可以计算出轨道参数,如离心率、倾角、长短轴等。

常用的轨道参数计算公式如下:1.离心率:e=√1+2El2μ(GM⊕)2)2.倾角:i=arccos(ℎzℎ3.长轴:a=−μT22E4.短轴:b=a√1−e2其中,E表示能量,l表示轨道角动量,μ表示标准引力参数,G表示引力常数,M⊕表示地球的质量,ℎ和ℎz分别表示轨道角动量和轨道角动量在z轴上的分量。

82第二节 龙格—库塔法

82第二节 龙格—库塔法
h k y x n O hk 1 k! 2
k
(1)
h h k 若令 yn1 y xn hy xn y xn y xn (2) 2! k! 则 y xn1 yn1 O hk 1
y0 k1 2 k2 hf x0 h 2, y0 k1 2
y0 k3
k4 hf x0 h, y0 k3
y0 k2 2 k3 hf x0 h 2, y0 k2 2
k
x1 x0 h y1 y0 k
数学学院 信息与计算科学系
0.1832292
0.1584376
数学学院 信息与计算科学系
接上图
0.4 0.5 0.5 0.6 0.6 1.341667 1.416026 1.412676 1.482627 1.483281 0.0745394 0.0710094 0.0708400 0.0673253
0.1416245
数学学院 信息与计算科学系
数学学院 信息与计算科学系
由表8-4可见,虽然四阶龙格-库塔方法每步要 计算四次 f 的值,但以h=0.2为步长ቤተ መጻሕፍቲ ባይዱ计算结果就
有5 位有效数字,而欧拉法与预估计-校正方法以
h=0.1为步长的计算结果才具有2 位与3 位有效数字.
如果步长 h 也取0.2,则结果的精度会更低.

即公式(2)为k 阶方法.
数学学院 信息与计算科学系
二、龙格-库塔方法(R-K方法)
R-K方法不是通过求导数的方法构造近似公式, 而是通过计算不同点上的函数值, 并对这些函数值作 线性组合, 构造近似公式, 再把近似公式与解的泰勒 展开式进行比较, 使前面的若干项相同 , 从而使近似 公式达到一定的阶数.

四阶龙格—库塔法的原理及其应用

四阶龙格—库塔法的原理及其应用

《四阶龙格—库塔法的原理及其应用》
龙格—库塔法(又称龙格库塔法)是由一系列有限的、独立的可能解组成的无穷序列,这些解中每个都与原来的数列相差一个常数。

它是20世纪30年代由匈牙利著名数学家龙格和库塔提出的,故得此名。

1.它的基本思想是:在n 阶方阵M 上定义一个函数,使得当n 趋于无穷时,它在m 中所表示的数值为M 的某种特征值,从而构造出一族具有某种特性的可计算函数f (x)= Mx+ C (其中C 为任意正整数)。

例如,若f (x)=(a-1) x+ C,则称之为(a-1) x 的龙格—库塔法。

2.它的应用很广泛,可以求解各类问题,且能将大量的未知数变换成少数几个已知数,因此它是近似计算的一种重要工具。

3.
它的优点主要有:(1)可以将多项式或不等式化成比较简单的形式;(2)对于同一问题可以用不同的方法来解决,并取得同样的结果;(3)适合处理高次多项式或者不等式,尤其适合处理多元函数的二次型。

第三部分龙格-库塔方法

第三部分龙格-库塔方法

内江师范学院数学与信息科学学院 吴开腾 制作
于是有
其中
y ( xn +1 ) − y ( xn ) = y '(ξ ), ξ ∈ ( xn , xn +1 ) h y ( xn +1 ) = y ( xn ) + hf (ξ , y (ξ ))
k * = f (ξ , y (ξ )) 称作区间 [ xn , xn +1 ] 上的平均斜率。 上的平均斜率 平均斜率。 问题:计算近似值y ( xn +1 ) 的关键是如何选择算法确定平均斜率 k *
(15)
f ( xn +1 , yn + h ( − k1 + 2 k 2 ))
内江师范学院数学与信息科学学院 吴开腾 制作
注释1 可以用Taylor展示证明格式(14) 注释1:可以用Taylor展示证明格式(14)具有三阶精 展示证明格式
度,并且还可以用类似的方法得到四阶及其以上的更高 阶精度的Runge-Kutta格式 阶精度的Runge-Kutta格式。 Runge 格式。
内江师范学院数学与信息科学学院 吴开腾 制作
h yn + ( k1 + 2 k 2 + 2 k3 +k 4) 6 f ( xn , y n ) h f ( x 1 , yn + k1 ) n+ 2 2 h f ( x 1 , yn + k 2 ) n+ 2 2 f ( xn +1 , yn + hk3 ) (16)
四阶龙格- 四阶龙格-库塔格式计算结果
xn
0.1 0.2 0.3 0.4 0.5
yn
欧拉格式计算结果 xn yn y ( xn )

矩阵特征值与特征向量的计算-Rung-Kutta方法

矩阵特征值与特征向量的计算-Rung-Kutta方法

每步须算Ki 的个数 2 3
4
5
6
可达到的最高精度 O(h2 ) O(h3 ) O(h4 ) O(h4 ) O(h5 )
7
O(h6 )
n8
O(hn−2 )
由于龙格-库塔法的导出基于泰勒展开,故精度主要受
解函数的光滑性影响。对于光滑性不太好的解,最好 采用低阶算法而将步长h 取小。
R − K方法的主要优缺点
二级R-K方法
二级R-K方法的形式为
其局部截断误差为 将 中的各项作Taylor展开
Taylor展开有 可得 令
二级R-K能达到的 最高阶数是二阶
常用的二级二阶R-K方法

,得
该方法称为改进的Euler公式(梯形公式的预估校正格式)

,得
该方法称为中点公式

,得
该方法称为Heun(休恩)方法
当 为实数时,得Euler法的绝对稳定区间是
二级R-K方法的绝对稳定区间
二阶二级R-K方法的计算公式为
由此可知,二阶二级R-K方法的绝对稳定区间是 当 为实数时,得绝对稳定区间是
一些常用方法的绝对稳定区间
R-K法的绝对稳定区域
k = 4 • 3. k =3
• 2. k=2 k = 1 • 1.
题相容的充分必要条件是该单步法至少是一阶方
法。
我们本章讨论的数值方法都是与原初值问题相容的!
收敛性
定义:对任意固定的 步法产生的解 ,均有
, 若初值问题的单
则称该方法是收敛的。
我们本章讨论的数值方法都是收敛的!
收敛性判别
定理7.3:设增量函数
在区域
上连续,并对变量y和h满足Lipschitz条件。如果单步 法与微分方程初值问题相容,则单步法收敛。

龙格库塔公式

龙格库塔公式

龙格-库塔(Runge-Kutta)公式是一种常用于数值解求微分方程的方法,其中包括了多种不同类型的龙格-库塔公式.
其中最常见的是四阶龙格-库塔公式(Runge-Kutta 4th order method),这种方法通过在当前时刻的函数值和当前时刻的导函数的近似值计算下一时刻的函数值。

它的计算公式如下:
k1 = hf(xn, yn)
k2 = hf(xn + h/2, yn + k1/2)
k3 = hf(xn + h/2, yn + k2/2)
k4 = hf(xn + h, yn + k3)
yn+1 = yn + (k1 + 2k2 + 2k3 + k4)/6
xn+1 = xn + h
其中k1,k2,k3,k4 是近似值,f(x,y)是微分方程,xn, yn 是当前时刻的函数值,xn+1, yn+1 是下一时刻的函数值, h是时间步长.
这种方法具有较高的精度,并且易于实现,但在求解某些类型的方程时,会存在稳定性问题.
除此之外,还有其他类型的龙格-库塔公式,比如二阶、三阶等等,它们也有各自的适用条件.。

龙格-库塔方法-文档资料

龙格-库塔方法-文档资料
3
c3a 2b32
c3a3 1
6
; 2
O (h4)
常见的2种三阶方法:
库塔三阶方法
h yn1yn6(k14k2k3)
k1
f(xn,yn);
k2
hh f(xn2,yn2k1)
k 3 f(x n h ,y n h k 1 2 h k 2 ) •5
四级方法:N = 4
局部截断误差 O ( h 5 )
可见误差随着 x n 的增加呈指数函数递减
当 f y 0 时,微分方程是不稳定的; 而 f y 0 时,微分方程是稳定的。
上面讨论的稳定性,与数值方法和方程中 f 有关
•21
实验方程: y y C ,R e () 0
D e f 3 对单步法 yn 1ynh(xn,yn,h )应用实验方程,
e n 1 e n h [ ( x n ,y ( x n ) , h ) ( x n ,y n , h ) ] T n 1
•15
因为单步法是 p 阶的:h0,0hh0满足|Tn1|Chp1
|e n 1| |e n| h L |e n| C h p 1|en |
其中 1hL,C hp1
•18
三、绝对稳定性 /*Absolute Stibility*/ 计算过程中产生的舍入误差对计算结果的影响
首先以Euler公式为例,来讨论一下舍入误差的传播:
yn1ynhf(xn,yn)
设实际计算得到的点 x n 的近似函数值为 yn yn n,
其中 y n 为精确值, n 为误差
yn1ynhf(xn,yn)
通过适当选取参数 1,2和 p 的值,使得公式具有 2阶精度!!
•3
由泰勒公式展开,要使公式具有 2 阶精度,只需

龙格库塔

龙格库塔

龙格-库塔(Runge-kutta)算法在MATLAB中的用法微分方程ode23:使用二阶龙格-库塔法求微分方程,调用格式为:[t,y]=ode23(‘fname’,tspan,y0)ode45:使用四阶龙格-库塔法求微分方程,调用格式为:[t,y]=ode45(‘fname’,tspan,y0)其中fname为由M函数定义的线性或者非线性微分方程的句柄函数名。

tspan形式为[t0,tf],表示求解区间。

y0是初始状态列向量。

t和y分别给出时间向量和相应的状态向量。

二阶龙格-库塔法(ode23):下面式2为Euler(欧拉法)增量,即一步起始端斜率,式3为一步终点端斜率。

所以式1仿真计算的增量实际上是取两点斜率的平均斜率来计算的,其精度高于Euler算法。

四阶龙格-库塔法(ode45):计算原理为预报-校正法,预报值采用Euler算出,下式又作了3次校正,因此计算精度远高于Euler算法。

e.g.二阶非线性系统的微分方程:x″ + 0.5*x′+ 2*x + x^2 = 0求系统在初始条件为x(0)=1,x′(0)=0的数值解。

建立M函数:function xdot=odetest(t,x)%龙格-库塔算法测试%二阶非线性系统微分方程xdotdot + 0.5*xdot + 2*x + x^2 = 0 %求系统在初始条件为x(0)=1,xdot(0)=0的数值解xdot=zeros(2,1);xdot(1)=x(2);xdot(2)=-0.5*x(2)-2*x(1)-x(1)^2;tspam=[0 20];x0=[0 1];[t,x]=ode23(@odetest,tspan,x0);plot(t,x)微分方程数值解曲线如下图:。

第二节_龙格_库塔方法

第二节_龙格_库塔方法

yn1 yn hf ( xn , yn )
n1 yn1 yn1
n1 n h[ f ( xn , yn ) f ( xn , yn )] [1 hf y ( xn , )]n
如果 |1 hf y | 1,则误差是不增的,故可认为是稳定的
例如:对于初值问题
y y
y(
x0
)
a
精确解为
k2
f ( xn
h 2
,
yn
h 2 k1)
hh k3 f ( xn 2 , yn 2 k2 )
k4 f ( xn h, yn hk3 )
例2:用经典的龙格-库塔方法
求解下列初值问题 h 0.1
dy 。 dx
y
2x y
x
(0,1)
解:经典的四阶龙格-库塔公式: y(0) 1
h
y
ae x x0
y y
而实际求解的初值问题为
y(
x0
)
a
a
精确解为 y (a a)e x x0在 xn 处的误差为ae xn x0
可见误差随着 xn的增加呈指数函数增长
y y
如果初值问题为
y( x0 ) a
精确解为 y ae x0 x
y y
实际求解的初值问题为
y( x0 )
2C( h)5 2
y( xn1)
y(h/ 2) n1
y( xn1)
y(h) n1
1 16

|
y(h/ 2) n1
y(h) n1
|
16( y( xn1)
y(h/ 2) n1
)
y( xn1)
y(h) n1
y( xn1)
y(h/ 2) n1

龙格库塔法-原理及程序实现

龙格库塔法-原理及程序实现

龙格库塔法一、基本原理:“龙格-库塔(Runge-Kutta)方法”是一种在工程上应用广泛的“高精度单步算法”。

由于此算法精度高,采取措施对误差进行抑制,可以得出四阶龙格-库塔公式,也就是在工程中应用广泛的经典龙格-库塔算法,即:yi+1=yi+h*( K1+ 2*K2 +2*K3+ K4)/6K1=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<stdio.h>#include<math.h>#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 20x0y0k1k2k3 k40.000000 2.000000-0.000000-0.500000-0.469238-0.8861310.250000 1.882308-0.885771-1.176945-1.129082-1.2800600.500000 1.599896-1.279834-1.295851-1.292250-1.2227280.750000 1.279948-1.228700-1.110102-1.139515-0.9901621.000000 1.000027-1.000054-0.861368-0.895837-0.7528521.2500000.780556-0.761584-0.645858-0.673410-0.5621891.5000000.615459-0.568185-0.481668-0.500993-0.4205371.7500000.492374-0.424257-0.361915-0.374868-0.3178552.0000000.400054-0.320087-0.275466-0.284067-0.2435982.2500000.329940-0.244935-0.212786-0.218538-0.1894822.5000000.275895-0.190295-0.166841-0.170744-0.1495632.7500000.233602-0.150068-0.132704-0.135399-0.1197033.0000000.200020-0.120024-0.106973-0.108868-0.0970483.2500000.172989-0.097256-0.087300-0.088657-0.0796183.5000000.150956-0.079757-0.072054-0.073042-0.0660303.7500000.132790-0.066124-0.060087-0.060818-0.0553054.0000000.117655-0.055371-0.050580-0.051129-0.0467434.2500000.104924-0.046789-0.042945-0.043363-0.0398334.5000000.094123-0.039866-0.036750-0.037072-0.0342024.7500000.084885-0.034226-0.031675-0.031926-0.029571xn=5.000000yn=0.076927。

数值分析9-3(龙格-库塔方法)

数值分析9-3(龙格-库塔方法)
语言实现龙格-库塔方法
总结词
除了Python和MATLAB,还有许多其他编 程语言可以用于实现龙格-库塔方法。
详细描述
例如C、Java和R等编程语言也提供了相应 的数值计算库或框架,可以实现龙格-库塔 方法。使用这些语言实现龙格-库塔方法需 要一定的编程基础和对相应语言的数值计算 库的了解。
龙格-库塔方法可以用于求解偏微分方程的数值解,通过将偏微分方程转化为常微分方程组,利用龙格 -库塔方法进行迭代求解,能够得到较为精确的结果。
积分方程的数值解
积分方程是描述函数与积分之间的关 系的数学模型,常见于物理、工程等 领域。
VS
龙格-库塔方法也可以用于求解积分 方程的数值解,通过将积分方程转化 为常微分方程组,利用龙格-库塔方 法进行迭代求解,能够得到较为精确 的结果。
重要性及应用领域
龙格-库塔方法是数值分析中非常重要的内容, 它为解决常微分方程提供了一种有效的数值方 法。
在科学、工程和经济学等领域中,许多问题都 可以转化为求解常微分方程的问题,因此龙格库塔方法具有广泛的应用价值。
例如,在物理学、化学、生物学、金融学等领 域中,龙格-库塔方法被广泛应用于模拟和预测 各种动态系统的行为。
数值分析9-3:龙格-库塔方法
目录
• 引言 • 龙格-库塔方法概述 • 龙格-库塔方法在数值分析中的应用 • 龙格-库塔方法的实现与编程 • 龙格-库塔方法的改进与优化 • 结论与展望
01 引言
主题简介
龙格-库塔方法是一种用于求解常微 分方程的数值方法。
它通过构造一个离散化的时间序列来 逼近微分方程的解,并利用已知的离 散点来计算新的离散点,逐步逼近微 分方程的真实解。
02 龙格-库塔方法概述
定义与原理

龙格-库塔方法(Runge-Kutta)

龙格-库塔方法(Runge-Kutta)

龙格-库塔方法〔Runge-Kutta〕3.2 Runge-Kutta法3.2.1 显式Runge-Kutta法的一般形式上节已给出与初值问题(1.2.1)等价的积分形式(3.2.1)只要对右端积分用不同的数值求积公式近似就可得到不同的求解初值问题(1.2.1)的数值方法,若用显式单步法(3.2.2)当,即数值求积用左矩形公式,它就是Euler法(3.1.2),方法只有一阶精度,若取(3.2.3)就是改进Euler法,这时数值求积公式是梯形公式的一种近似,计算时要用二个右端函数f的值,但方法是二阶精度的.若要得到更高阶的公式,则求积分时必须用更多的f值,根据数值积分公式,可将(3.2.1)右端积分表示为注意,右端f中还不能直接得到,需要像改进Euler法(3.1.11)一样,用前面已算得的f值表示为(3.2.3),一般情况可将(3.2.2)的 表示为(3.2.4)其中这里均为待定常数,公式(3.2.2),(3.2.4)称为r级的显式Runge-Kutta法,简称R-K方法.它每步计算r个f值(即由前面(i-1)个已算出的表示,故公式是显式的.例),而ki如当r=2时,公式可表示为(3.2.5)其中.改进Euler 法(3.1.11)就是一个二级显式R-K 方法.参数取不同的值,可得到不同公式.3.2.2 二、三级显式R-K 方法对r=2的显式R-K 方法(3.2.5),要求选择参数,使公式的精度阶p 尽量高,由局部截断误差定义11122211()()[(,())(,)]n n n n n n n T y x y x h c f x y x c f x a h y b hk ++=--+++ (3.2.6)令,对(3.2.6)式在处按Taylor 公式展开,由于将上述结果代入(3.2.6)得要使公式(3.2.5)具有的阶p=2,即,必须(3.2.7)即由此三式求的解不唯一.因r=2,由〔3.2.5〕式可知,于是有解(3.2.8)它表明使(3.2.5)具有二阶的方法很多,只要都可得到二阶精度R-K方法.若取,则,则得改进Euler法(3.1.11),若取,则得,此时(3.2.5)为(3.2.9)其中称为中点公式.改进的Euler法(3.1.11)与中点公式(3.2.9)是两个常用的二级R-K方法,注意二级R-K方法只能达到二阶,而不可能达到三阶.因为r=2只有4个参数,要达到p=3则在(3.2.6)的展开式中要增加3项,即增加三个方程,加上(3.2.7)的三个方程,共计六个方程求4个待定参数,验证得出是无解的.当然r=2,p=2的R-K方法(3.2.5)当取其他数时,也可得到其他公式,但系数较复杂,一般不再给出.对r=3的情形,要计算三个k值,即其中将按二元函数在处按Taylor公式展开,然后代入局部截断误差表达式,可得可得三阶方法,其系数共有8个,所应满足的方程为(3.2.10)这是8个未知数6个方程的方程组,解也是不唯一的,通常.一种常见的三级三阶R-K方法是下面的三级Kutta方法:(3.2.11)附:R-K 的三级Kutta 方法程序如下function y = DELGKT3_kuta(f, h,a,b,y0,varvec) format long; N = (b-a)/h;y = zeros(N+1,1); y(1) = y0; x = a:h:b;var = findsym(f); for i=2:N+1K1 = Funval(f,varvec,[x(i-1) y(i-1)]);K2 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]); K3 = Funval(f,varvec,[x(i-1)+h y(i-1)-h*K1+K2*2*h]);y(i) = y(i-1)+h*(K1+4*K2+K3)/6; %满足c1+c2+c3=1,(1/6 4/6 1/6)endformat short; 3.2.3 四阶R-K 方法与步长的自动选择利用二元函数Taylor 展开式可以确定(3.2.4)中r=4,p=4的R-K 方法,其迭代公式为111223344()n n y y h c k c k c k c k +=++++其中1(,)n n k f x y =,2221(,(,))n n n n k f x a h y b hf x y =++,而33311322(,)n n k f x a h y b hk b hk =+++ 44411422433(,)n n k f x a h y b hk b hk b hk =++++共计13个参数待定,Taylor 展开分析局部截断误差,使得精度达到四阶,即误差为5()O h 。

龙格库塔法的基本原理

龙格库塔法的基本原理

龙格库塔法的基本原理嘿,朋友们!今天咱来唠唠龙格库塔法的基本原理。

咱就说啊,这龙格库塔法就像是一个超级厉害的解题高手!它能帮我们在面对那些复杂得让人头疼的数学问题时找到答案。

你想想,有时候那些数学式子就像一团乱麻,让你无从下手。

但龙格库塔法呢,它就有办法一点点地把这团乱麻给解开。

它通过巧妙地计算和推测,一步一步地接近问题的真相。

就好比你要去一个陌生的地方,你不知道该怎么走。

龙格库塔法就像是一个聪明的向导,它能给你指出一条清晰的路。

它不是那种随随便便指个方向就完事儿的,而是经过深思熟虑,考虑了各种因素之后才做出的决定。

龙格库塔法在很多领域都大显身手呢!比如在物理、工程这些地方,那可真是立下了汗马功劳。

没有它,好多难题都没法解决,好多伟大的工程可能都没法实现。

它的原理其实也不复杂,就是通过一系列的计算步骤来逼近真实的答案。

但可别小看这些步骤哦,它们就像一个个小小的拼图,组合起来就能呈现出一幅完整的画面。

比如说,它会在不同的点上进行计算,就像在不同的地方观察一样,然后把这些观察到的信息综合起来,得出一个比较准确的结果。

这是不是很神奇?而且啊,龙格库塔法还特别灵活。

它可以根据不同的问题进行调整和优化,就像一个武林高手能根据对手的情况随时改变招式一样。

你说,这么厉害的龙格库塔法,咱能不好好了解了解吗?咱得知道它是怎么工作的,怎么帮我们解决问题的呀!不然岂不是浪费了这么好的一个工具?总之,龙格库塔法就是数学世界里的一颗璀璨明星,照亮了我们解决难题的道路。

它让那些看似不可能的问题变得有可能,让我们在知识的海洋里能更加自由地航行。

所以啊,大家可别小瞧了它哟!。

滤波算法 龙格库塔算法-概述说明以及解释

滤波算法 龙格库塔算法-概述说明以及解释

滤波算法龙格库塔算法-概述说明以及解释1.引言1.1 概述概述:滤波算法和龙格库塔算法是计算机科学领域中常用的算法之一,它们在数据处理和数值计算中有着重要的应用价值。

滤波算法被广泛应用于信号处理、图像处理、通信系统等领域,用于消除信号中的噪声和提高数据的质量。

而龙格库塔算法则是一种常用的数值求解微分方程的方法,能够有效地对复杂的数学模型进行数值求解,具有较高的准确性和稳定性。

本文将分别介绍滤波算法和龙格库塔算法的原理、优缺点以及应用领域,希望读者通过本文能够对这两种算法有更深入的了解,并在实际应用中能够灵活运用。

1.2 文章结构本文将分为四个部分来探讨滤波算法和龙格库塔算法。

首先在引言部分,对滤波算法和龙格库塔算法进行简要介绍,并说明本文的结构和目的。

接着在第二部分,详细介绍滤波算法的概念、常见算法和应用场景,以便读者对滤波算法有个全面的了解。

然后在第三部分,深入探讨龙格库塔算法的简介、原理和优缺点,帮助读者更好地理解这一种数值计算方法。

最后,在结论部分对两种算法进行总结,并展望未来可能的发展方向,以及得出结论。

通过以上四个部分的内容,读者能够全面了解和掌握滤波算法和龙格库塔算法的相关知识。

1.3 目的本文的主要目的是介绍和探讨滤波算法和龙格库塔算法这两种在计算机科学和工程领域中广泛应用的算法。

通过对这两种算法的概述、原理和应用进行详细分析,能够帮助读者全面了解它们的工作原理和特点。

同时,通过对这两种算法的比较和讨论,可以帮助读者更好地理解它们在不同应用场景下的适用性和优缺点。

此外,本文还旨在为读者提供一个深入学习和掌握这两种算法的基础知识和入门指南。

通过本文的学习,读者可以加深对滤波算法和龙格库塔算法的理解,为进一步的研究和实践打下坚实的基础。

同时,希望本文能够激发读者对算法领域的兴趣,促使他们深入研究和探索更多先进的算法及其应用。

2.滤波算法2.1 滤波算法概述滤波算法是一种用于处理信号或数据的技术,其主要目的是通过去除噪声或不需要的信息,从而提取出所需的信号或数据。

龙格库塔

龙格库塔

数值分析中,龙格-库塔法(Runge-Kutta)是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。

这些技术由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。

经典四阶龙格库塔法龙格库塔法的家族中的一个成员如此常用,以至于经常被称为“RK4”或者就是“龙格库塔法”。

令初值问题表述如下。

则,对于该问题的RK4由如下方程给出:其中这样,下一个值(y n+1)由现在的值(y n)加上时间间隔(h)和一个估算的斜率的乘积决定。

该斜率是以下斜率的加权平均:∙k1是时间段开始时的斜率;∙k2是时间段中点的斜率,通过欧拉法采用斜率k1来决定y在点t n + h/2的值;∙k3也是中点的斜率,但是这次采用斜率k2决定y值;∙k4是时间段终点的斜率,其y值用k3决定。

当四个斜率取平均时,中点的斜率有更大的权值:RK4法是四阶方法,也就是说每步的误差是h5阶,而总积累误差为h4阶。

注意上述公式对于标量或者向量函数(y可以是向量)都适用。

显式龙格库塔法显示龙格-库塔法是上述RK4法的一个推广。

它由下式给出其中如果要求方法有精度p则还有相应的条件,也就是要求舍入误差为O(h p+1)时的条件。

这些可以从舍入误差本身的定义中导出。

例如,一个2阶精度的2段方法要求b1 + b2 = 1, b2c2 = 1/2, 以及b2a21 = 1/2。

在Matlab下输入:edit,然后将下面两行百分号之间的内容,复制进去,保存%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%function dxdt=ode_Miss_ghost(t,x)%分别用x(1),x(2),x(3),x(4)代替N1,P1,N2,P2N1=x(1);P1=x(2);N2=x(3);P2=x(4);K=2;tau_c=3e-9;tan_p=6e-12;beta =5e-5;delta=0.692;eta =0.0001;fm =8e6;Ith =26e-3;Ib =1.5*Ith;Im =0.3*Ith;I1=Ib+Im*sin(2*pi*fm*t)+K*P2;I2=Ib+Im*sin(2*pi*fm*t)+K*P1;dxdt=[(I1/Ith-N1-(N1-delta)/(1-delta)*P1)/tau_e;((N1-delta)/(1-delta)*(1-eta*P1)*P1-P1+beta*N1)/tau_p;(I2/Ith-N2-(N2-delta)/(1-delta)*P2)/tau_e;((N2-delta)/(1-delta)*(1-eta*P2)*P2-P2+beta*N2)/tau_p;]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%在Matlab下面输入:t_start=0;t_end=2e-9;y0=[1e-3;1e-4;0;0]; %初值[x,y]=ode15s('ode_Miss_ghost',[0,t_end],y0);plot(x,y);legend('N1','P1','N2','P2');xlabel('x');Mathlab定步长龙格-库塔法[345阶]、自适应步长rkf45经常看到很多朋友问定步长的龙格库塔法设置问题,下面吧定步长三阶、四阶、五阶龙格库塔程序贴出来,有需要的可以看看ODE3 三阶龙格-库塔法CODE:function Y = ode3(odefun,tspan,y0,varargin)%ODE3 Solve differential equations with a non-adaptive method of order 3.% Y = ODE3(ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN] integrates% the system of differential equations y' = f(t,y) by stepping from T0 to% T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector.% The vector Y0 is the initial conditions at T0. Each row in the solution% array Y corresponds to a time specified in TSPAN.%% Y = ODE3(ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters% P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...).%% This is a non-adaptive solver. The step sequence is determined by TSPAN% but the derivative function ODEFUN is evaluated multiple times per step.% The solver implements the Bogacki-Shampine Runge-Kutta method of order 3.%% Example% tspan = 0:0.1:20;% y = ode3(@vdp1,tspan,[2 0]);% plot(tspan,y(:,1));% solves the system y' = vdp1(t,y) with a constant step size of 0.1, % and plots the first component of the solution.%if ~isnumeric(tspan)error('TSPAN should be a vector of integration steps.');endif ~isnumeric(y0)error('Y0 should be a vector of initial conditions.');endh = diff(tspan);if any(sign(h(1))*h <= 0)error('Entries of TSPAN are not in order.')endtryf0 = feval(odefun,tspan(1),y0,varargin{:});catchmsg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr];error(msg);endy0 = y0(:); % Make a column vector.if ~isequal(size(y0),size(f0))error('Inconsistent sizes of Y0 and f(t0,y0).');endneq = length(y0);N = length(tspan);Y = zeros(neq,N);F = zeros(neq,3);Y(:,1) = y0;for i = 2:Nti = tspan(i-1);hi = h(i-1);yi = Y(:,i-1);F(:,1) = feval(odefun,ti,yi,varargin{:});F(:,2) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,1),varargin{:});F(:,3) = feval(odefun,ti+0.75*hi,yi+0.75*hi*F(:,2),varargin{:});Y(:,i) = yi + (hi/9)*(2*F(:,1) + 3*F(:,2) + 4*F(:,3));endY = Y.';ODE4 四阶龙格-库塔法CODE:function Y = ode4(odefun,tspan,y0,varargin)%ODE4 Solve differential equations with a non-adaptive method of order 4.% Y = ODE4(ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN] integrates % the system of differential equations y' = f(t,y) by stepping from T0 to% T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector.% The vector Y0 is the initial conditions at T0. Each row in the solution% array Y corresponds to a time specified in TSPAN.%% Y = ODE4(ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters % P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...).%% This is a non-adaptive solver. The step sequence is determined by TSPAN % but the derivative function ODEFUN is evaluated multiple times per step.% The solver implements the classical Runge-Kutta method of order 4.%% Example% tspan = 0:0.1:20;% y = ode4(@vdp1,tspan,[2 0]);% plot(tspan,y(:,1));% solves the system y' = vdp1(t,y) with a constant step size of 0.1,% and plots the first component of the solution.%if ~isnumeric(tspan)error('TSPAN should be a vector of integration steps.');endif ~isnumeric(y0)error('Y0 should be a vector of initial conditions.');endh = diff(tspan);if any(sign(h(1))*h <= 0)error('Entries of TSPAN are not in order.')endtryf0 = feval(odefun,tspan(1),y0,varargin{:});catchmsg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr]; error(msg);endy0 = y0(:); % Make a column vector.if ~isequal(size(y0),size(f0))error('Inconsistent sizes of Y0 and f(t0,y0).');endneq = length(y0);N = length(tspan);Y = zeros(neq,N);F = zeros(neq,4);Y(:,1) = y0;for i = 2:Nti = tspan(i-1);hi = h(i-1);yi = Y(:,i-1);F(:,1) = feval(odefun,ti,yi,varargin{:});F(:,2) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,1),varargin{:}); F(:,3) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,2),varargin{:}); F(:,4) = feval(odefun,tspan(i),yi+hi*F(:,3),varargin{:});Y(:,i) = yi + (hi/6)*(F(:,1) + 2*F(:,2) + 2*F(:,3) + F(:,4));endY = Y.';定步长RK4(自编):CODE:function Y=RungeKutta4(f,xn,y0)% xn=0:.1:1;% y0=1;% y_n=[];% f=@(X1,Y1) Y1-2*X1/Y1;y_n=[];h=diff(xn(1:2));for i=1:length(xn)-1K1=f(xn(i),y0);K2=f(xn(i)+h/2,y0+h*K1/2);K3=f(xn(i)+h/2,y0+h*K2/2);K4=f(xn(i)+h,y0+h*K3);y_n=[y_n;y0+h/6*(K1+2*K2+2*K3+K4)];y0=y_n(end);endY=y_n;ODE5 五阶龙格-库塔法CODE:function Y = ode5(odefun,tspan,y0,varargin)%ODE5 Solve differential equations with a non-adaptive method of order 5.% Y = ODE5(ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN] integrates % the system of differential equations y' = f(t,y) by stepping from T0 to% T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector.% The vector Y0 is the initial conditions at T0. Each row in the solution% array Y corresponds to a time specified in TSPAN.%% Y = ODE5(ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters % P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...).%% This is a non-adaptive solver. The step sequence is determined by TSPAN % but the derivative function ODEFUN is evaluated multiple times per step.% The solver implements the Dormand-Prince method of order 5 in a general % framework of explicit Runge-Kutta methods.%% Example% tspan = 0:0.1:20;% y = ode5(@vdp1,tspan,[2 0]);% plot(tspan,y(:,1));% solves the system y' = vdp1(t,y) with a constant step size of 0.1,% and plots the first component of the solution.if ~isnumeric(tspan)error('TSPAN should be a vector of integration steps.');endif ~isnumeric(y0)error('Y0 should be a vector of initial conditions.');endh = diff(tspan);if any(sign(h(1))*h <= 0)error('Entries of TSPAN are not in order.')endtryf0 = feval(odefun,tspan(1),y0,varargin{:});catchmsg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr];error(msg);endy0 = y0(:); % Make a column vector.if ~isequal(size(y0),size(f0))error('Inconsistent sizes of Y0 and f(t0,y0).');endneq = length(y0);N = length(tspan);Y = zeros(neq,N);% Method coefficients -- Butcher's tableau%% C | A% --+---% | BC = [1/5; 3/10; 4/5; 8/9; 1];A = [ 1/5, 0, 0, 0, 03/40, 9/40, 0, 0, 044/45 -56/15, 32/9, 0, 019372/6561, -25360/2187, 64448/6561, -212/729, 0 9017/3168, -355/33, 46732/5247, 49/176, -5103/18656];B = [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84]; % More convenient storageA = A.';B = B(:);nstages = length(B);F = zeros(neq,nstages);Y(:,1) = y0;for i = 2:Nti = tspan(i-1);hi = h(i-1);yi = Y(:,i-1);% General explicit Runge-Kutta frameworkF(:,1) = feval(odefun,ti,yi,varargin{:});for stage = 2:nstageststage = ti + C(stage-1)*hi;ystage = yi + F(:,1:stage-1)*(hi*A(1:stage-1,stage-1));F(:,stage) = feval(odefun,tstage,ystage,varargin{:});endY(:,i) = yi + F*(hi*B);endY = Y.';-------------------------------------------------------------------------------------------------------------------自适应步长RKF45(相当于ode45)ODE45 是4阶方法提供候选解,5阶方法控制误差。

矩阵龙格库塔

矩阵龙格库塔

要求电流就是求解矩阵微分方程:(R+pM(t))*I(t)+M(t)*pI(t)-U(t)=0,其中p是求导,R是6*6常数矩阵,M(t)是6*6的时变矩阵,U(t)是6*1的时变矩阵,求I(t),也是6*1的矩阵。

已知条件:M=[0,0,0, -15727/10000*sin(5/12*pi+80*pi*t)*pi,15727/10000*cos(1/4*pi+80*pi*t)*pi, 15 727/10000*sin(1/12*pi+80*pi*t)*pi;0,0,0,15727/10000*sin(1/12*pi+80*pi*t)*pi, -15727/100 00*sin(5/12*pi+80*pi*t)*pi,15727/10000*cos(1/4*pi+80*pi*t)*pi;0,0,0,15727/10000*cos(1/4 *pi+80*pi*t)*pi,15727/10000*sin(1/12*pi+80*pi*t)*pi, -15727/10000*sin(5/12*pi+80*pi*t)*p i;-15727/10000*sin(5/12*pi+80*pi*t), 15727/10000*sin(1/12*pi+80*pi*t)*pi, 15727/100 00*cos(1/4*pi+80*pi*t)*pi, 0,0, 0;15727/10000*cos(1/4*pi+80*pi*t)*pi, -15727/10000*sin(5/12*pi+80*pi*t)*pi, 15727/1 0000*sin(1/12*pi+80*pi*t)*pi, 0,0, 0;15727/10000*sin(1/12*pi+80*pi*t)*pi, 15727/10000*cos(1/4*pi+80*pi*t)*pi, -15727/ 10000*sin(5/12*pi+80*pi*t)*pi, 0,0, 0];U=[380*cos(100*pi*t);380*cos(100*pi*t-2*pi/3);380*cos(100*pi*t+2*pi/3);0;0;0];R=[0.0247,0,0,0,0,0;0,0.0247,0,0,0,0;0,0,0.0247,0,0,0;0,0,0,0.0193,-0.0193,0;0,0,0,0,0.0193,-0.0193;0,0,0,1,1,1];电流I的初值是:I(0)=[0;0;0;0;0;0];t是时间预期最后的结果是6个电流与时间的关系谢谢分享到:举报| 2010-06-22 16:47提问者采纳global R M Usyms tR=[0.0247,0,0,0,0,0;0,0.0247,0,0,0,0;0,0,0.0247,0,0,0;0,0,0,0.0193,-0.0193,0;0,0,0,0,0.0193,-0.0193;0,0,0,1,1,1];M=[0,0,0,-15727/10000*sin(5/12*pi+80*pi*t)*pi,15727/10000*cos(1/4*pi+80*pi*t)*pi,15727/10 000*sin(1/12*pi+80*pi*t)*pi;0,0,0,15727/10000*sin(1/12*pi+80*pi*t)*pi, -15727/10000*sin(5/12*pi+80*pi*t)*pi, 15727/10000*cos(1/4*pi+80*pi*t)*pi;0,0,0, 15727/10000*cos(1/4*pi+80*pi*t)*pi,15727/10000*sin(1/12*pi+80*pi*t)*pi, -15727/10000*sin(5/12*pi+80*pi*t)*pi;-15727/10000*sin(5/12*pi+80*pi*t),15727/10000*sin(1/12*pi+80*pi*t)*pi,15727/10000*cos(1/4*pi+80*pi*t)*pi,0,0,0;15727/10000*cos(1/4*pi+80*pi*t)*pi,-15727/10000*sin(5/12*pi+80*pi*t)*pi,15727/10000*sin(1/12*pi+80*pi*t)*pi,0,0,0;15727/10000*sin(1/12*pi+80*pi*t)*pi, 15727/10000*cos(1/4*pi+80*pi*t)*pi, -15727/10000*sin(5/12*pi+80*pi*t)*pi,0,0,0];U=[380*cos(100*pi*t);380*cos(100*pi*t-2*pi/3);380*cos(100*pi*t+2*pi/3);0;0;];DM=diff(M,t); %就是M对t求导%%%要调节的参数在这里%%注意,你的M有点奇异,计算很快发散掉了,你检察一下相关的参数吧。

龙格库塔积分算法

龙格库塔积分算法

龙格库塔法龙格库塔法是常用于模拟常微分方程的解的重要的一类隐式或显式迭代法。

这些技术由数学家C. Runge和M.W. Kutta于1900年左右发明。

由于此算法精度高,采取措施对误差进行抑制,所以其实现原理也较复杂。

龙格库塔法是一种在工程上应用广泛的高精度单步算法,可以应用在物理、工程、控制、动力学中,如模糊控制、弹道分析以及分析光纤特性等,在系统仿真中得到广泛应用。

龙格库塔法源自于相应的泰勒级数方法,在每一插值节点用泰勒级数展开,其截断误差阶数也是,根据可省略更高阶的导数计算, 这种方法可构造任意阶数的龙格库塔法。

其中4 阶龙格库塔法是最常用的一种方法。

因为它相当精确、稳定、容易编程。

在计算中一般不必使用高阶方法, 因为附加的计算误差可由增加精度来弥补。

如果需要较高的精度, 可采取减小步长的方法即可。

4 阶龙格库塔法的精度类似4 阶泰勒级数法的精度。

1、初值问题对于一阶常微分方程的初值问题根据常微分方程的理论可知,此初值问题的解在区间[a,b]上存在,且唯一。

2、离散化取步长h=(b-a)/n,将区间[a , b]分成n个子区间:a=<=b在其中任意两点的曲线段上,根据积分中值定理,一段光滑曲线上至少有一点,它的斜率与整段曲线的平均斜率相同,得=y’() (0<<1)其中,=可以将上式改写成y()=y()+h*K (2.1)其中K为平均斜率,K=f()公式(2.1)表明,如果能够确定平均斜率K,就可以根据(2.1)式得到y()的值。

欧拉法和龙格库塔法就是用不同方法确定不同精度的平均斜率K,从而求得y()的近似值。

3、Euler法欧拉法虽然精度低,但它是最简单的一种显式单步法,也是龙格库塔法的基础。

首先,令、为y() 及y()的近似值,并且令平均斜率K=f(),即以点的斜率作为平均斜率K,便得到欧拉公式=+h* f() (3.1)4、改进的欧拉法此种方法是取、两点的斜率的平均值作为平均斜率K,即K= ,其中、均为y()以及y()的近似值,就得到改进后的欧拉公式(4.1)其中、分别为、两点的斜率值,即= ,=在上面的(4.1)式中,k2是未知的,采用一种叫预报法的方法来求解。

用龙格_库塔法求解非线性方程组

用龙格_库塔法求解非线性方程组

目录用龙格- 库塔法求解非线性方程组 (1)1引言 (1)2 龙格- 库塔法的基本原理 (1)3 多启动延拓法 (2)3 微分同伦 (3)4 显式龙格- 库塔(Runge - Kutta) 法 (3)参考文献 (4)用龙格- 库塔法求解非线性方程组摘 要 本文介绍了一种求解非线性方程组的新方法龙格- 库塔法。

关键词 多启动延拓法 非线性方程组 龙格- 库塔法 1引言如何求解非线性方程组: )1()...,,(,0)(21nT n R x x x x x F ∈==在非线性科学技术中有重要的作用,这里T n f f f F )...,,(21=是从n 维欧氏空间n R 到n R 的光滑映射。

大家知道,Newton 迭代法和它的变体是解非线性方程组的经典方法,至今都是基本的重要的。

但是它的半径一般很小,需要很好的初值0x 。

因此,Newton 法仅当初值0x 很接近真解x '时才是收敛的。

但由于真解x '事先不知道,实际上寻找一个很好的初值0x 本身就是很困难的。

其重要的研究成果见文献[ 1 ]和[ 2 ] 。

因此在20 世纪70 年代,发展了一种延拓法(ex2tension method ,EM) (也称同伦算法) ,它有大范围的收敛性,取初值的容许范围被显著扩大了。

于是,它成为解非线性方程组最有效的方法,其重要的研究成果见文献[ 1 ] 。

但将EM 解某些复杂问题时,仍常常发散,因EM 仍严重依赖初值0x ,使延拓曲线严重弯曲所致。

在某种程度上为了克服这些缺点,本文改用多启动延拓法(multi - starting extend method ,MSEM) ,寻找很接近真值的初值,然后,采用微分同伦利用龙格- 库塔法(Runge - Kutta ,method) 求出非线性方程组的解x '。

2 龙格- 库塔法的基本原理假设函数Y 是微分方程),(Y X f Y ='满足初始条件00)(Y X Y =的解。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

要求电流就是求解矩阵微分方程:(R+pM(t))*I(t)+M(t)*pI(t)-U(t)=0,
其中p是求导,
R是6*6常数矩阵,
M(t)是6*6的时变矩阵,
U(t)是6*1的时变矩阵,
求I(t),也是6*1的矩阵。

已知条件:
M=[0,0,0, -15727/10000*sin(5/12*pi+80*pi*t)*pi,15727/10000*cos(1/4*pi+80*pi*t)*pi, 15 727/10000*sin(1/12*pi+80*pi*t)*pi;0,0,0,15727/10000*sin(1/12*pi+80*pi*t)*pi, -15727/100 00*sin(5/12*pi+80*pi*t)*pi,15727/10000*cos(1/4*pi+80*pi*t)*pi;0,0,0,15727/10000*cos(1/4 *pi+80*pi*t)*pi,15727/10000*sin(1/12*pi+80*pi*t)*pi, -15727/10000*sin(5/12*pi+80*pi*t)*p i;-15727/10000*sin(5/12*pi+80*pi*t), 15727/10000*sin(1/12*pi+80*pi*t)*pi, 15727/100 00*cos(1/4*pi+80*pi*t)*pi, 0,
0, 0;
15727/10000*cos(1/4*pi+80*pi*t)*pi, -15727/10000*sin(5/12*pi+80*pi*t)*pi, 15727/1 0000*sin(1/12*pi+80*pi*t)*pi, 0,
0, 0;
15727/10000*sin(1/12*pi+80*pi*t)*pi, 15727/10000*cos(1/4*pi+80*pi*t)*pi, -15727/ 10000*sin(5/12*pi+80*pi*t)*pi, 0,
0, 0];
U=[380*cos(100*pi*t);
380*cos(100*pi*t-2*pi/3);
380*cos(100*pi*t+2*pi/3);
0;
0;
0];
R=[0.0247,0,0,0,0,0;
0,0.0247,0,0,0,0;
0,0,0.0247,0,0,0;
0,0,0,0.0193,-0.0193,0;
0,0,0,0,0.0193,-0.0193;
0,0,0,1,1,1];
电流I的初值是:I(0)=[0;0;0;0;0;0];
t是时间
预期最后的结果是6个电流与时间的关系
谢谢
分享到:
举报| 2010-06-22 16:47提问者采纳
global R M U
syms t
R=[
0.0247,0,0,0,0,0;
0,0.0247,0,0,0,0;
0,0,0.0247,0,0,0;
0,0,0,0.0193,-0.0193,0;
0,0,0,0,0.0193,-0.0193;
0,0,0,1,1,1
];
M=[
0,0,0,-15727/10000*sin(5/12*pi+80*pi*t)*pi,15727/10000*cos(1/4*pi+80*pi*t)*pi,15727/10 000*sin(1/12*pi+80*pi*t)*pi;
0,0,0,15727/10000*sin(1/12*pi+80*pi*t)*pi, -15727/10000*sin(5/12*pi+80*pi*t)*pi, 15727/10000*cos(1/4*pi+80*pi*t)*pi;
0,0,0, 15727/10000*cos(1/4*pi+80*pi*t)*pi,15727/10000*sin(1/12*pi+80*pi*t)*pi, -15727/10000*sin(5/12*pi+80*pi*t)*pi;
-15727/10000*sin(5/12*pi+80*pi*t),15727/10000*sin(1/12*pi+80*pi*t)*pi,
15727/10000*cos(1/4*pi+80*pi*t)*pi,0,0,0;
15727/10000*cos(1/4*pi+80*pi*t)*pi,
-15727/10000*sin(5/12*pi+80*pi*t)*pi,15727/10000*sin(1/12*pi+80*pi*t)*pi,0,0,0;
15727/10000*sin(1/12*pi+80*pi*t)*pi, 15727/10000*cos(1/4*pi+80*pi*t)*pi, -15727/10000*sin(5/12*pi+80*pi*t)*pi,0,0,0
];
U=[
380*cos(100*pi*t);
380*cos(100*pi*t-2*pi/3);
380*cos(100*pi*t+2*pi/3);
0;
0;
];
DM=diff(M,t); %就是M对t求导
%%%要调节的参数在这里
%%注意,你的M有点奇异,计算很快发散掉了,你检察一下相关的参数吧。

%%det(subs(M,t,0))
I0=[0;0;0;0;0;0];
tstart=0;
tend=0.1;
dt=0.1;
%%%end
tout=tstart:dt:tend;
n=length(tout);
M_t_dt=subs(M,t,tstart); %P=subs(P,'t',x)就是把P表达式中所有't',都用具体的x值代替U_t_dt=subs(U,t,tstart);
DM_t_dt=subs(DM,t,tstart);
II=I0;
for i=1:n-1
tt=tout(i);
M_t=M_t_dt;
U_t=U_t_dt;
DM_t=DM_t_dt;
M_t_dt_2=subs(M,t,tt+dt/2);
U_t_dt_2=subs(U,t,tt+dt/2);
DM_t_dt_2=subs(DM,t,tt+dt/2);
M_t_dt=subs(M,t,tt+dt);
U_t_dt=subs(U,t,tt+dt);
DM_t_dt=subs(DM,t,tt+dt); %(R+pM(t))*I(t)+M(t)*pI(t)-U(t)=0
k1=dt*M_t \(U_t -(R+DM_t )*(II(:,end) ));
k2=dt*M_t_dt_2\(U_t_dt_2-(R+DM_t_dt_2)*(II(:,end)+0.5*k1));
k3=dt*M_t_dt_2\(U_t_dt_2-(R+DM_t_dt_2)*(II(:,end)+0.5*k2));
k4=dt*M_t_dt \(U_t_dt -(R+DM_t_dt )*(II(:,end)+k3 ));
I_t=II(:,end)+(k1+2*k2+2*k3+k4)/6; II=[II,I_t];
end
plot(tout',II')。

相关文档
最新文档