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

算法原理和程序框图

龙格—库塔法是一种求其准确解y(x)在一系列点xi处y(xi)的近似值yi的方法,yi称为数值解。经典的四阶龙格库塔法方程如下:

y'=f(t,y),y(t0)=y0输出按如下求解yn+1=yn+h(k1+2k2+3k3+4k4)/6其中                       k1=f(tn,yn)

k2=f(tn+h/2,yn+hk1/2)

k3=f(tn+h/2,yn+hk2/2)

k4=f(tn+h,yn+hk3)

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



  • 0
    点赞
  • 30
    收藏
    觉得还不错? 一键收藏
  • 1
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值