Matlab 龙格库塔解常微分方程组练习

练习。使用四阶龙格库塔法求解常微分方程组，通用性较佳。附加一个振动方程求解的案例。振动方程是一个二阶微分方程，转化为两个方程组以后用编写的代码求解。

%自动变步长四阶龙格库塔法求常微分方程的解，微分方程的格式应该是
%y‘（x）=f（xy）; 初始值为（x0y0）
%返回值为x、y两个向量
function [xy]=runge4_adaptStep（fx0xNy0）
    x（1）=x0;
    y（1）=y0;
    h=（xN-x0）/100;%假定一个初始步长，在计算中根据情况，这个值可能会被程序更改
    i=1;
    x_temp=x（1）;
    
    stepHalfed=0;%记录h是否被减半过。对于某些情况，或许会出现以h为步长时步长太大，但以h/2为步长时步长又太小
                    %这种情况下，为了避免步长在h和h/2之间来回跳跃进入死循环，将采用h/2为步长
    e2=1e-6;
    e1=1e-5;
    while x（i）        y1=getY（fx_tempy（i）h）;%以h为步长计算，得到x_temp+h处的近似值
        y11=getY（fx_tempy（i）h/2）;%以h/2为步长连续计算两次，得到x_temp+h处的近似值
        y12=getY（fx_temp+h/2y11h/2）;
        e=abs（（y1-y12）/y1）;%计算上面两种途径求得的x+h处的y值之间的比值
        if e>e1%步长太大，减小步长后重新计算一轮
            h=h/2;
            stepHalfed=1;
        elseif e>e2%步长合适，保存得到的值并准备下一点的计算
            x（i+1）=x_temp+h;
            y（i+1）=y1;
            x_temp=x（i+1）;
            i=i+1;
            stepHalfed=0;
        else%步长太小
            if stepHalfed==1%此时步长是经过缩小后的，如果步长放大，又会过大，会陷入死循环。
                                %因此接受当前步长，保存当前步结果且准备进入下一点的计算
                x（i+1）=x_temp+h;
                y（i+1）=y1;
                x_temp=x（i+1）;
                i=i+1;
                stepHalfed=0;
            else%此时步长不是经过缩小后的，确实是太小，需要放大
                h=2*h;%二话不说，把步长放大，从头再试试                
            end
        end        
    end
end

function [newY]=getY（fxyh）
    k1=f（xy）;
    k2=f（x+0.5*hy+0.5*h*k1）;
    k3=f（x+0.5*hy+0.5*h*k2）;
    k4=f（x+hy+h*k3）;
    newY=y+h/6*（k1+2*k2+2*k3+k4）;
end

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        1833  2018-12-27 00:19  runge4_adaptStep.m
     文件         492  2018-12-27 00:19  test_equationSet.m