CP080-计算物理常微分方程解资料

k

hf
( xn

1 2
h,
yn

2k3 1 2 k)

k4
)
公式推导 见P37-44

k3

hf
( xn

1 2
h,
yn

1 2
k2 )
k4 hf (xn h, yn k3 )
例：用R-K方法解例题8.1.1
例..
dQ Q dt RC
Q() Q
f
(
xn1
,
y(k) n
)]
（k ,,,..）（. 3）
当|
y (k ) n

yn(k )
|

时, 取yn

y (k ) n
8.2 改进Euler方法
对于(2)计算yn1 ,由于迭代工作量较大, 一般只 迭代一次,构成一类预估 校正算法,即

y(0) n1
%1、预报 Q(i)=Q(i-1)+h*I(i-1); I(i)=I(i-1)+h*f(Q(i-1),I(i-1)); %2、计算Q k1=h*I(i-1); k2=h*I(i); Q(i)=Q(i-1)+1/2*(k1+k2); %3、计算I k1=h*f(Q(i-1),I(i-1)); k2=h*f(Q(i),I(i)); I(i)=I(i-1)+1/2*(k1+k2); end plot(t,Q,'r',t,I,'m');%改进欧拉法的曲线
hf
(tn

h, Qn

k3 )

h(
Qn k RC
)
function main %龙格库塔：一阶常微分方程 global E R C; E=10;R=10;C=0.01;h=0.01;t=[0:h:1];Q(1)=E/C; for i=2:length(t)
k1=h*f(Q(i-1)); k2=h*f(Q(i-1)+k1/2); k3=h*f(Q(i-1)+k2/2); k4=h*f(Q(i-1)+k3); Q(i)=Q(i-1)+1/6*(k1+2*k2+2*k3+k4); end plot(t,Q,‘b’);hold on%龙格库塔方法的曲线 plot(t,Q(1)*exp(-1*t./(R*C)),'.');%理论曲线
function y=f(Q,I) global E R C L; y=(E-Q/C-I*R)/L;
（用改进的Euler法解）：
、先用Euler法预报：
Q n

Qn

hIn
,
(n

,,,...)
I n

In

h(E
Qn
/
C

I n R)
/
L, (n

,,,...)
、改进的Euler法：
function y=f(x) global E R C; y=-1*x./(R*C);
例：用R-K方法解例题8.1.2
例..2

dQ dt

I
L dI IR Q E
dt
C
Q() I()


Qn

k

k
hf hf
dy dx
|xn
y( xn
h) h
y(xn )

f (xn , y(xn ))
以yn表示y(xn )的近似值，则有
yn1 yn hf (xn , yn ) （n 0,1,2,..）. (1)
这就是显式的Euler公式，它可以从y0出发，逐次 算出y1, y2 , y3...。
不相等，但一般取成相等的，这时 h b a 。在这些节点上采用离散化 n
方法（通常用泰勒展开），将上述初值问题化成关于离散变量的相应问题。 把这个相应问题的解 yn 作为 y(xn)的近似值。这样求得的 yn 就是上述初值 问题在节点 xn 上的数值解。一般说来，不同的离散化导致不同的方法。

Qn


k hf
Qn

1 6
(k

2k

2k3
(t
n
,
Qn
)

h(
Qn RC
)

k4
)
k

hf
(tn

1 2
h, Qn

1 2
k )

h(
Qn
k RC
/
)

k3

hf
(tn

1 2
h, Qn

1 2
k2
)

h(
Qn
k RC
/
)
k4


Qn

Qn


(k

k )
k


k


h h
dQn dt dQn dt
hIn h(
I
n
)

(n

,,,...)


I
n

In


(k

k )
和k


k
h dIn
dt

h
dI
n
dt
h（E h（E
常微分方程的数值解法
绪论
在工程和科学计算中，所建立的各 种常微分方程的初值或边值问题，除很 少几类的特殊方程能给出解析解，绝大 多数的方程是很难甚至不可能给出解析 解的，其主要原因在于积分工具的局限 性。因此，人们转向用数值方法去解常 微分方程，并获得相当大的成功，讨论 和研究常微分方程的数值解法是有重要 意义的。
常微分方程描写的物理现象
镭的衰变规律 单摆的运动 RLC振荡电路 物理场计算


l
d
dt
m
d
dt

mg sin
dQ dt

I
L
dI dt
IR

Q C

E
常微分方程
我们考虑一阶常微分方程初值问题
dy dx

f (x,
y)
y(x ) y
QE(i)=QE(i-1)+h*f(QE(i-1));%欧拉法 k1=h*f(QEG(i-1));%改进欧拉法 k2=h*f(QEG(i-1)+k1); QEG(i)=QEG(i-1)+1/2*(k1+k2); end plot(t,QE,'r');%欧拉法的曲线 hold on plot(t,QEG,'b');%改进欧拉法的曲线 plot(t,Q0*exp(-1*t./(R*C)),‘.’);%理论曲线 legend('欧拉法的曲线','改进欧拉法的曲线', '理论曲线')
y(x)

y(x0 )
y'(x0 )(x
x0 )
y"(x0 ) 2!
(x
x0 )2

y(n) (x0 ) n!
(x
x0 )n
...
xx0xxn8n .1h Euler方法
如果用xn1代替x0 ,于是该式可离散为：
y( xn
h) h
y(xn )

f
( xn1 ,
y(x)

y(x0
)

y' ( x0 )(x

x0 )

y" ( x0 ) 2!
(x

x0 )2

y(n) (x0 ) n!
(x

x0 )n
...
x x8n .1h Euler方法
x0 xn
为实现这一目标，Euler方法首先将微分算子离
散化，并用xn代替x0 , 于是该式可离散为：
Qn (tn (tn
,I
1 6 (k 2k
n ) hIn
1
1
2 h, In 2
2k3 k4 ) l) h(In
1 2
l )
k3

hf
(tn

1 2
h,
In

1 2
l2 )

h(In

1 2
l
)
k4

hf
(tn

h, In
l3 )

h(In
Qn / C Qn
In R）/ L / C InR）/
L

(n

,,,...)
function main %改进欧拉法：二阶常微分方程 global E R C L; E=10;R=100;C=0.01;L=10; Q0=0;I0=0;h=0.01;t=[0:h:10];Q(1)=Q0;I(1)=I0; for i=2:length(t)
function y=f(x) global E R C; y=-1*x./(R*C);
例8.1.2如图，当开关 K接通时，试分别用 Euler 法和改进 Euler 法求解求 RLC 电路中C的电量和电流的变化情 况.
K
R=100欧
E=10V
L=10H
C=10mF
解：分析电路的常微分方程组：
dQ

Qn1
改进的Euler法：k1

hf
Qn (Qn )
1 2
(k1 k2 ) h( Qn )
RC

k2

hf
(Qn

k1 )

h(
Qn k1 ) RC
(n 0,1,2,...)
function main global E R C; E=10;R=10;C=0.01; Q0=E/C; h=0.01; t=[0:h:1]; QE(1)=Q0; QEG(1)=Q0; for i=2:length(t)
y ( xn1 ))
以yn表示y(xn )的近似值，则有
yn1 yn hf (xn1, yn1) （n 0,1,2,..）. (2)
这就是隐式的Euler公式或向后Euler方法，它与显式
的不同在于，它每算一步要解函数方程(2)才能得到
yn
。
1
8.1 Euler方法-梯形公式
function y=f(Q,I) global E R C L; y=(E-Q/C-I*R)/L;
8.3 龙格-库塔（R-K）方法
思想：取多点处斜率的加权平均为平均斜 率，从而减小误差。
四阶公式：

1 yn yn 6 (k 2k
k hf (xn , yn )
在区间[a, b]上的解，其中 f (x, y)为 x, y 的已知 函数，y0 为给定的初始值，将上述问题的精确 解记为 y(x)。
常微分方程数值解基本思想
数值方法的基本思想是：在解的存在区间上取 n + 1 个节点 a x0 x1 x2 xn b
这里差 hi xi1 xi ，i = 0,1, …, n 称为由 xi 到 xi+1 的步长。这些 hi 可以

1 2
l )
dQ

dt

I
dI (E - IR - Q ) / L
dt
C
Q()
如果取以上两式的算术平均值的结果，则得
yn

yn

h[f
(xn , yn )
f
(xn1, yn1)]
（n ,,,..）.
ห้องสมุดไป่ตู้称为梯形公式。
计算yn时常用以下迭代式：

y(k) n

yn
hf (xn , yn )

y (k ) n

yn

h[f
(xn , yn )
2K 1
R=10欧
C=0.01法
E=10伏
解：分析得电量满足如下常微分方程

dQ dt


Q RC
Q(t0 ) Q0
初始状态：
Q0 E / C 10 / 0.01; 取h 0.01, 计算t [0,1]上结果，此时
Euler法：Qn1

Qn

h(
Qn ) RC
(n 0,1,2,...)

yn

hf
(xn ,
yn )

y (1) n1`

yn

h[ 2
f
(xn ,
yn )
f
( xn1 ,
y(0) n1
)]
并取yn1 yn(1)1。
8.2 改进Euler方法
上式还常写成

yn k
yn hf (xn ,

(k
yn )

k
)

dt

I
L dI IR Q E
dt
C
Q(0) 0
I(0) 0

初始状态：Q0 0; I0 0; 取h 0.01, 计算t [0,10]上结果，此时 Euler法：
Qn1 Qn hIn , (n 0,1,2,...)
In1 In h(E Qn / C InR) / L, (n 0,1,2,...)
function main %欧拉法：二阶常微分方程 global E R C L; E=10;R=100;C=0.01;L=10; Q0=0;I0=0;h=0.01;t=[0:h:10]; Q(1)=Q0; I(1)=I0; for i=2:length(t)
Q(i)=Q(i-1)+h*I(i-1);%Q I(i)=I(i-1)+h*f(Q(i-1),I(i-1));%I end plot(t,Q,'r',t,I,'m');%欧拉法的曲线
k hf (xn h, yn k)
(n ,,,...)
该式称为改进Euler方法, 亦可写成
yn

yn

h[
f
( xn

yn )
f
( xn ,
yn
hf
(xn ,
yn ))]
例8.1.1 如图，当开关从1接至2时，求RC电路 放电过程中电量的变化情况, 试分别用Euler法 和改进Euler法求解。