有限差分法

基本概念

总的思想就是把微分方程变成差分方程（离散化），用精度换取计算的可进行性。

• 几种差分方法及其误差：
  (来源：https://wenku.baidu.com/view/aa580fd619e8b8f67d1cb919.html）
  
  
  
  
  weixin_43793265 提醒上图左边应该是 f(x+h) + f(x-h)

MATLAB实现

以一个对时间进行微分的微分方程组为例（数模国赛2014A）：

其中
$$\dot{m}=\frac{F}{\nu }$$
各参数取值（都使用国际单位制，故单位不再给出）：
$$\begin{gathered} G=6.672\times 10^{-11}; \\ M=7.3477\times 10^{22}; \\ F=7500; \\ \nu =2940; \\ {m}_0=2.4\times 10^3; \end{gathered}$$
已知条件：
$$\begin{gathered} v_{x0}=1692.46; \\ {v}_{y0}=0; \\ {r}_0=1749.37\times 10^3; \end{gathered}$$
使用前向差分可把方程化为：
$$\{\begin{array}{c}{v}_{xn+1}=(a_{n}cos{\beta }_n-\frac{v_{xn}{v}_{yn}}{r_n})h+v_{xn}\\ v_{yn+1}=(-\frac{GM}{r_n^2}+\frac{v_{xn}^2}{r_n}+a_{n}sin{\beta }_n)h+v_{yn}\\ r_{n+1}=h{y}_{yn}+r_n\\ a_n=\frac{F}{m_0-\frac{F}{\nu }nh}\\ tan{\beta }_n=\frac{v_{yn}}{v_{xn}}\end{array}$$
用以下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