计算机应与用基础-5-常微分方程求解 .ppt

龙格-库塔法是求解常微分方程较常用的一种方法，它通过巧妙的线性组合，在显式格式的情况下获得理想的计算精度，大大提高了计算速度。该方法的推导过程不是本课程要研究的对象，本课程主要研究其实际应用,，故直接给出各类龙格-库塔公式。 1、 二阶龙格-库塔 3、四阶龙格—库塔公式 龙格-库塔方法——步长的选择 怎样选取合适的步长，这在实际计算中是很重要的。由于步长愈小，每步计算的截断误差就愈小；但在一定的求解范围内，需要完成的步数就愈多，这不但引起计算量的增加，而且计算步数的增加往往造成舍入误差的严重积累，反而降低了计算精度。上面介绍的龙格-库塔方法是对定步长（即步长h为常数）而言的，但为了保证精确度，一种有效的措施是在计算过程中自动进行步长调整，即变步长技巧。 把它与上式相比，可得 经整理可得 这表明以 作为y(xn+1) 的近似值，其误差可用先后两次计算结果之差来表示，因而，只需考察 是否成立。若成立，则可将 作为y(xn+1)的近似值；若不成立，则将步长再次对折进行计算，直到不等式成立为止，并取最后的 作为计算结果。以上方法就是计算过程中自动选择步长的方法，也称为变步长方法。从表面上看，为了选择适当的步长，每一步的计算量增加了，但从总体考虑，工作量往往还是减少的。 下面用一个例子说明： 用标准四阶龙格-库塔方法计算得 由积分法算得y(2)=2.23607,当h=0.5时相差0.00013，而h=0.25误差小于0.000001，但当 h=0.125时虽然足够精确但计算次数却比h=0.25多了一倍。因此合理选择步长既能保证精度又能减少计算量。 5.5 高阶常微分方程求解 Function [x,y]=fdiff(funcs, tspan, x0f,n) t0=tspan(1); tfinal=tspan(2); ga=x0f(1); gb=x0f(2); for i=1:n x(i)=t0+h*(i-1); end x0=x(1:n-1); t=-2+h^2*feval(funcs,x0,2); tmp= feval(funcs,x0,1); v=1+h*tem/2; w=1-h*tem/2; b= h^2*feval(funcs,x0,3); b(1)=b(1)-w(1)*ga; b(n-1)=b(n-1)-v(n-1)*gb; b=b’; A=diag(t); for i=1:n-2 A(i,i+1)=v(i); A(i+1,i)=w(i+1); end y=inv(A)*b; x=[x, tfinal]; y=[ga; y; gb]’; 5.4 微分方程的解析解方法 边值问题求解之三-符号算法 y=dsolve(f1,f2,…,fm) y=dsolve(f1,f2,…,fm,’x’) % 指明自变量 求 x’(t)=x(t)(1-x2(t))的解析解 Syms x t X=dsolve(Dx=x*(1-x^2)) 例5-8 X = [ 1/(1+exp(-2*t)*C1)^(1/2)] [ -1/(1+exp(-2*t)*C1)^(1/2)] 【例5-5】输入信号 syms t y u=exp(-5*t)*cos(2*t+1)+5; uu=5*diff(u,t,2)+4*diff(u,t)+2*u? uu =87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10 5.4 微分方程的解析解方法 求以下微分方程的通解： 步骤1 步骤2 Syms t y y=dsolve([‘D4y+10*D3y+35*D2y+50*Dy+24*y=’, ... ‘87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10’]) y = -87/52*cos(1)*cos(t)^2*exp(-5*t)-87/260*cos(1)*sin(t)*cos(t)*exp(-5*t)+87/104*exp(-5*t)*cos(1)+23/26*exp(-5*t)*sin(1)-87/520*sin(1)*cos(2