不能无限制地减小h的值的两条原因: (1)减慢计算速度 (2)增加累积误差 在对微分方程求解过程中应采取的三个措施 (1)选择适当的步长 (2)改进近似算法精度 (3)采用变步长方法
7.2.2四阶定步长 Runge-Kuta算法 及 MATLAB实现 四阶定步长 Runge-Kutta算法: k1=hf(tk, k k2=f( k +-wk xI k3=hf tk+3,k+ k4=hf(k+h,k+k3)其中h为计算步长 下一个步长的状态变量值为: xk+1=k+(k1+2k2+2k3+k4)
7.2.2 四阶定步长Runge-Kutta算法 及 MATLAB 实现
function[ tout, yout-rk4( defile, tspan,y0)%y0初值列向量 to=tspan(1); th=tspan(2) if length (tspan)<=3, h=tspan ( 3); tspan=to, th, h] else,h- tspan(2)- tspan(1);th+ tspan(end);end%等间距数组 tout=[to: h: th, yout=L for t=tout k1=h*eval(lodefile(t,yo defile是一个字符串变 量,为表示微分方程的文件名 k2=h*eval(lodefile(t+h/,y0+0.5*k1)D k3=h*eval(lodefile '(t+h/2, y0+0.5*k2) D k4=h*eval(lodefile (t+h, y0+k3)] y0=y0+(k1+2*k2+2*k3+k4)/6; yout=Lout; yoT end %由效果看,该算法不是一个较好的方法
function [tout,yout]=rk_4(odefile,tspan,y0) ％y0初值列向量 t0=tspan(1); th=tspan(2); if length(tspan)<=3, h=tspan(3); % tspan=[t0,th,h] else, h=tspan(2)-tspan(1); th=tspan(end); end ％等间距数组 tout=[t0:h:th]'; yout=[]; for t=tout' k1=h*eval([odefile ‘(t,y0)’]); % odefile是一个字符串变 量，为表示微分方程的文件名。 k2=h*eval([odefile '(t+h/2,y0+0.5*k1)']); k3=h*eval([odefile '(t+h/2,y0+0.5*k2)']); k4=h*eval([odefile '(t+h,y0+k3)']); y0=y0+(k1+2*k2+2*k3+k4)/6; yout=[yout; y0']; end ％由效果看，该算法不是一个较好的方法
7.2.3一阶微分方程组的数值解 7.2.3.1四阶五级 Runge-Kutta-Felhberg算法 假设当前的步长为kk,定义6个k变量: ki- hkf tk+a hk, k+>Biik 1,2,…,6 下一步的状态向量:ck+1=k+yk 定义一个误差向量 k (1-y;)k 通过误差向量调节步长,此为自动变步长方法 四阶五级RKF算法有参量系数表
7.2.3 一阶微分方程组的数值解 7.2.3.1 四阶五级Runge-Kutta-Felhberg算法 通过误差向量调节步长，此为自动变步长方法。 四阶五级RKF算法有参量系数表
7.2.3.2基于 MATLAB的微分方程 求解函数 格式1:直接求解 Lt, x]=ode45(Fun, [to, t1, xo) 格式2:带有控制参数 Lt, x]=ode45(Fun, [to, tl, Xo, options 格式3:带有附加参数 Lt, x]=ode45(Fun, [to, tl, Xo, options, p1, p2, . . [to,td求解区间,x初值问题的初始状态变量
7.2.3.2 基于 MATLAB 的微分方程 求解函数 格式1： 直接求解 [t,x]=ode45(Fun,[t0,tf],x0) 格式2： 带有控制参数 [t,x]=ode45(Fun,[t0,tf],x0，options) 格式3： 带有附加参数 [t,x]=ode45(Fun,[t0,tf],x0,options,p1,p2,…) [t0,tf]求解区间， x0初值问题的初始状态变量