四阶龙格-库塔法求解常微分方程的初值问题

算法原理和程序框图

龙格—库塔法是一种求其准确解y(x)在一系列点x_i处y（x_i）的近似值y_i的方法，y_i称为数值解。经典的四阶龙格库塔法方程如下：

y'=f(t,y),y(t₀)=y₀输出按如下求解y_n+1=y_n+h(k₁+2k₂+3k₃+4k₄)/6其中                       k₁=f(t_n,y_n)

k₂=f(t_n+h/2,y_n+hk₁/2)

k3=f(t_n+h/2,y_n+hk₂/2)

k4=f(t_n+h,y_n+hk₃)

4.2程序使用说明

本程序单独编写了一阶初值问题单独求解和高阶常微分方程的求解问题。常微分方程的阶数，x的上下限和步长直接输入就好。常微分方程则在系统提示input后再输入，并且按照y(1,1)=y,y（2,1）=y’,y(3,1)=y’’的格式输入,然后输入初始条件f(x),f’(x),f’’(x)等计算结果。最后求得所要的函数值。

4.3计算实例展示

例题为课本304页计算实习9.2。计算y（1）的结果如下：

ans =

  -1.000000000000000

  -0.847436458333333

  -0.689482936121419

  -0.525724908067489

  -0.355719511320969

  -0.178993729355764

   0.004957535888930

   0.196673510288970

   0.396729719312157

   0.605740292781225

   0.824360410550626

   1.053288897488131

   1.293270976642268

   1.545101189994835

   1.809626496745704

   2.087749559656611

% ** 文件名:Runge_Kutta.m
% 
% ** 创建人:**
% 
% ** 日 期:2016.12.14
% 
% ** 描 述:四级四阶龙格－库塔法求解常微分方程的初值问题
% 
% ** 函数：第304 页计算实习 9.2   
% 
% ** 参考教材：《数值分析》李乃成，梅立泉，科学出版社
%
clear
clc
format long
% % % m=input('请输入常微分方程的阶数m=');
% % % a=input('请输入x下限a=');
% % % b=input('请输入x上限b=');
% % % h=input('请输入步长h=');
% % % ym=input('令y(1,1)=y,y（2,1）=y’,y(3,1)=y’’...请输入ym=','s');     %输入的时候必须按照这个形式输入y1=y(1,1);

m=3;
a=0;
b=1;
h=0.05;
%ym=y(3,1)+y(2,1)-y(1,1)+2*x-3;
ym=input('令y(1,1)=y,y（2,1）=y’,y(3,1)=y’’...请输入ym=','s');     %输入的时候必须按照这个形式输入y1=y(1,1);
% if m==1                                               %一阶初值问题单独求解
%     mm=(b-a)/h;
%     y(1,1)=input('请输入在初值点的函数值f(a)=');
%     x=a;
%     y11(1)=y(1,1);
%     for k1=2:(mm+1)
%         y1=y(1,1);
%         K(1,1)=h*(eval(ym));                         %计算K1
%         x=x+h/2;
%         y(1,1)=y1+K(1,1)/2;
%         y1=y(1,1);
%         K(1,2)=h*(eval(ym));                         %计算K2
%         x=x;
%         y(1,1)=y1+K(1,2)/2-K(1,1)/2;
%         y1=y(1,1);
%         K(1,3)=h*(eval(ym));                          %计算K3
%         x=x+h/2;
%         y(1,1)=y1+K(1,3)-K(1,2)/2;
%         y1=y(1,1);
%         K(1,4)=h*(eval(ym));                          %计算K4
%         y11(k1)=y11(k1-1)+(K(1,1)+2*K(1,2)+2*K(1,3)+K(1,4))/6;
%         y(1,1)=y11(k1);
%         x=a+(k1-1)*h;
%        
%     end
% y11
% else                                                 %高阶初值问题
    mm=(b-a)/h;                                     %一共要求解mm个数据点
    for k2=1:m                                      %读取初值条件
        fprintf('请输入%d阶导数的初值f(%d)(a)=\n',(k2-1),(k2-1));
        y(k2,1)=input('=');
    end
    for k2=1:m                                        
         y22(1,k2)=y(k2,1);                          %先把初值保存在矩阵y22(m,n)中,m表示第几个所求点，n表示第n阶初值
    end
%%%%初值

% y22(1,1)=-1;
% y22(1,2)=3;
% y22(1,3)=2;
    x=a;
    for k4=2:(mm+1)                                   %求解mm个数据点的循环
        for k=1:(m-1)                                 %计算K1，包括每一阶的K1
            K(k,1)=h*y(k+1,1);                        %y(k+1,1)中k+1表示第k+1阶，1表示第一个点；K(k,1)中k表示阶数，1表示K1
        end
        K(m,1)=h*(eval(ym));
        x=x+h/2;                                      %求解K1之前，先重新对x和y赋值
        for k3=1:m              
            y(k3,1)=y(k3,1)+K(k3,1)/2;
        end
        for k=1:(m-1)                                  %计算K2
            K(k,2)=h*y(k+1,1);
        end
        K(m,2)=h*(eval(ym));
        x=x;
        for k3=1:m
            y(k3,1)=y(k3,1)-K(k3,1)/2+K(k3,2)/2;
        end
        for k=1:(m-1)                                  %计算K3
            K(k,3)=h*y(k+1,1);
        end
        K(m,3)=h*(eval(ym));
        x=x+h/2;
        for k3=1:m
            y(k3,1)=y(k3,1)+K(k3,3)-K(k3,2)/2;  %这里容易出错
        end
        for k=1:(m-1)                                     %计算K4
            K(k,4)=h*y(k+1,1);
        end
        K(m,4)=h*(eval(ym));
        for k5=1:m
            y22(k4,k5)=y22(k4-1,k5)+(K(k5,1)+2*K(k5,2)+2*K(k5,3)+K(k5,4))/6;         %这里，除了要求出下一个点的数值，还要求出相应的导数值
        end
        for k6=1:m                                       %除了对y(1,1)重新赋值外，还要对y(2,1)等重新赋值
            y(k6,1)=y22(k4,k6); 
        end                                                      
        x=a+(k4-1)*h;
    end
        y22(:,1)
%end