基本概念
总的思想就是把微分方程变成差分方程(离散化),用精度换取计算的可进行性。
- 几种差分方法及其误差:
(来源:https://wenku.baidu.com/view/aa580fd619e8b8f67d1cb919.html)
weixin_43793265 提醒上图左边应该是 f(x+h) + f(x-h)
MATLAB实现
以一个对时间进行微分的微分方程组为例(数模国赛2014A):
其中
m
˙
=
F
ν
\dot{m}=\frac{F}{\nu}
m˙=νF
各参数取值(都使用国际单位制,故单位不再给出):
G
=
6.672
×
1
0
−
11
;
M
=
7.3477
×
1
0
22
;
F
=
7500
;
ν
=
2940
;
m
0
=
2.4
×
1
0
3
;
G=6.672\times 10^{-11};\\ M=7.3477\times 10^{22};\\ F=7500;\\ \nu =2940;\\ m_0=2.4\times 10^3;
G=6.672×10−11;M=7.3477×1022;F=7500;ν=2940;m0=2.4×103;
已知条件:
v
x
0
=
1692.46
;
v
y
0
=
0
;
r
0
=
1749.37
×
1
0
3
;
v_{x0}=1692.46;\\ v_{y0}=0;\\ r_0=1749.37\times 10^3;
vx0=1692.46;vy0=0;r0=1749.37×103;
使用前向差分可把方程化为:
{
v
x
n
+
1
=
(
a
n
c
o
s
β
n
−
v
x
n
v
y
n
r
n
)
h
+
v
x
n
v
y
n
+
1
=
(
−
G
M
r
n
2
+
v
x
n
2
r
n
+
a
n
s
i
n
β
n
)
h
+
v
y
n
r
n
+
1
=
h
y
y
n
+
r
n
a
n
=
F
m
0
−
F
ν
n
h
t
a
n
β
n
=
v
y
n
v
x
n
\left \{\begin{matrix} v_{xn+1}=(a_n cos\beta_n-\dfrac{v_{xn}v_{yn}}{r_n})h+v_{xn}\\ v_{yn+1}=(-\dfrac{GM}{r_n^2}+\dfrac{v_{xn}^2}{r_n}+a_nsin\beta_n)h+v_{yn}\\ r_{n+1}=hy_{yn}+r_n\\ a_n=\dfrac{F}{m_0-\frac{F}{\nu}nh}\\ tan\beta_n=\dfrac{v_{yn}}{v_{xn}} \end{matrix}\right.
⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧vxn+1=(ancosβn−rnvxnvyn)h+vxnvyn+1=(−rn2GM+rnvxn2+ansinβn)h+vynrn+1=hyyn+rnan=m0−νFnhFtanβn=vxnvyn
用以下matlab程序计算
%% main.m
clear;
h=0.1; %步长
N=100000; %最大迭代次数
%初始化
vy=zeros(1,N);
vx=zeros(1,N);
r=zeros(1,N);
%初始值
vy(1)=0;
vx(1)=1692.46;
r(1)=1749.37*10^3;
%开始递推
for n=1:N-1
[vy(n+1),vx(n+1),r(n+1)]=DiTuiFun(vy(n),vx(n),r(n),h,n);
if((r(n+1)-1737013)<=12000 || vx(n+1)<=0)
end_num=n+1;
break;
end
end
r_reslut=r(1:end_num);
vx_reslut=vx(1:end_num);
vy_reslut=vy(1:end_num);
figure(1);
plot(1:1:end_num,r_reslut);
figure(2);
plot(1:1:end_num,vx_reslut);
figure(3);
plot(1:1:end_num,vy_reslut);
%% DiTuiFun.m
function [vynn,vxnn,rnn]=DiTuiFun(vyn,vxn,rn,h,n)
%参数
G=6.672*10^-11;
M=7.3477*10^22;
F=7500;
miu=2940;
m0=2.4*10^3;
%计算
v=sqrt(vyn^2+vxn^2);
sin_beta=vyn/v;
cos_beta=-vxn/v;
an=F/(m0-F*n*h/miu);
rnn=h*vyn+rn;
vxnn=(an*cos_beta-vxn*vyn/rn)*h+vxn;
vynn=(-G*M/rn^2+vxn^2/rn+an*sin_beta)*h+vyn;
end