当前位置: 首页 > 计算模拟 >微分方程组的四阶龙格库塔公式求解matlab版

微分方程组的四阶龙格库塔公式求解matlab版

作者 longwen8
来源: 小木虫 350 7 举报帖子
+关注

微分方程组的四阶龙格库塔公式,求解matlab版的代码。下面是我的代码,运行不出来。
function varargout=rungekutta(varargin)
clc,clear
x0=0;xn=1.0;y0=1;h=0.05;%h为步长
[y,x]=rungekutta4(x0,xn,y0,h);
Function z=f(x,y);
z=y-2*x/y;
function [y,x]=rungekutta4(x0,xn,y0,h)
x=x0:h:xn;
n=(xn-x0)/h;
y1=x;
y1(1)=y0;
for i=1:n %龙格库塔法的算法
    K1=f(x(i),y1(i));
    K2=f(x(i)+h/2,y1(i)+h/2*K1);
    K3= f(x(i)+h/2,y1(i)+h/2*K2);
K4= f(x(i)+h,y1(i)+h*K3);
y1(i+1)=y1(i)+h/6*(K1+2*K2+2*K3+K4);
    yy1(i)=(1+2*x(i+1)) ^0.5;
    error(i)=y1(i+1)-yy1(i);
    yy(i+1)=yy1(i);
    e(i+1)=error(i);
end
   y=y1;
n=(xn-x0)/h;
fprintf(‘i       x(i)     yi数值解     y(i)真值    误差\n’);
fprintf(‘-------------------------------------------------------------\n’);
for i=1:n
    fprintf(‘%2d %12.4f %14.8f %14.8f %12.8f\n’,i, x(i+1), y(i+1), yy(i+1), error(i));
    end

    plot(x,e,‘r.’);%画出误差曲线
    xlable(‘x轴’),ylable(‘误差’);
    title(‘步长为0.05时的误差曲线’); 返回小木虫查看更多

今日热帖
  • 精华评论
  • hzlhm

    一、几个自定义函数位置错误。应为这样放置
    1、主程序 ,
    function rungekutta()
    clc,clear
    x0=0;xn=1.0;y0=1;h=0.05;%h为步长
    [y,x]=rungekutta4(x0,xn,y0,h);
    end
    2、 四阶龙格库塔函数程序
    function [y,x]=rungekutta4(x0,xn,y0,h)
    。。。
    end
    3、z函数程序
    function z=f(x,y);
    z=y-2*x/y;
    end
    二、xlable(‘x轴’),ylable(‘误差’);这句代码中函数书写错误。应为
    xlabel('x轴'),ylabel('误差');
    三、修改后运行得到如下结果

  • longwen8

    引用回帖:
    2楼: Originally posted by hzlhm at 2020-12-18 18:12:48
    一、几个自定义函数位置错误。应为这样放置
    1、主程序 ,
    function rungekutta()
    clc,clear
    x0=0;xn=1.0;y0=1;h=0.05;%h为步长
    =rungekutta4(x0,xn,y0,h);
    end
    2、 四阶龙格库塔函数程序
    function =rungeku ...

    什么结果?大神。

  • hzlhm

    引用回帖:
    3楼: Originally posted by longwen8 at 2020-12-18 19:43:08
    什么结果?大神。...

    运行结果:
    i       x(i)       yi数值解        y(i)真值      误差
    -------------------------------------------------------------
    1       0.0500     1.04880886     1.04880885   0.00000001
    2       0.1000     1.09544514     1.09544512   0.00000002
    3       0.1500     1.14017546     1.14017543   0.00000004
    4       0.2000     1.18321600     1.18321596   0.00000005
    5       0.2500     1.22474493     1.22474487   0.00000006
    6       0.3000     1.26491113     1.26491106   0.00000007
    7       0.3500     1.30384056     1.30384048   0.00000008
    8       0.4000     1.34164088     1.34164079   0.00000009
    9       0.4500     1.37840498     1.37840488   0.00000011
    10       0.5000     1.41421368     1.41421356   0.00000012
    11       0.5500     1.44913781     1.44913767   0.00000014
    12       0.6000     1.48323985     1.48323970   0.00000015
    13       0.6500     1.51657526     1.51657509   0.00000017
    14       0.7000     1.54919353     1.54919334   0.00000019
    15       0.7500     1.58113904     1.58113883   0.00000021
    16       0.8000     1.61245178     1.61245155   0.00000023
    17       0.8500     1.64316793     1.64316767   0.00000026
    18       0.9000     1.67332034     1.67332005   0.00000028
    19       0.9500     1.70293895     1.70293864   0.00000031
    20       1.0000     1.73205115     1.73205081   0.00000034
    图像上传不上来

  • longwen8

    引用回帖:
    2楼: Originally posted by hzlhm at 2020-12-18 18:12:48
    一、几个自定义函数位置错误。应为这样放置
    1、主程序 ,
    function rungekutta()
    clc,clear
    x0=0;xn=1.0;y0=1;h=0.05;%h为步长
    =rungekutta4(x0,xn,y0,h);
    end
    2、 四阶龙格库塔函数程序
    function =rungeku ...

    function rungekutta()
    clc,clear
    x0=0;xn=1.0;y0=1;h=0.05;%h为步长
    [y,x]=rungekutta4(x0,xn,y0,h);
    End
    function [y,x]=rungekutta4(x0,xn,y0,h)
    x=x0:h:xn;
    n=(xn-x0)/h;
    y1=x;
    y1(1)=y0;
    for i=1:n %龙格库塔法的算法
        K1=f(x(i),y1(i));
        K2=f(x(i)+h/2,y1(i)+h/2*K1);
        K3= f(x(i)+h/2,y1(i)+h/2*K2);
    K4= f(x(i)+h,y1(i)+h*K3);
    y1(i+1)=y1(i)+h/6*(K1+2*K2+2*K3+K4);
        yy1(i)=(1+2*x(i+1))^0.5;
        error(i)=y1(i+1)-yy1(i);
        yy(i+1)=yy1(i);
        e(i+1)=error(i);
    end
    function z=f(x,y);
    z=y-2*x/y;
    end
       y=y1;
    n=(xn-x0)/h;
    fprintf('i       x(i)     yi数值解     y(i)真值    误差\n');
    fprintf('-------------------------------------------------------------\n');
    for i=1:n
        fprintf(‘%2d %12.4f %14.8f %14.8f %12.8f\n’,i, x(i+1), y(i+1), yy(i+1), error(i));
        end
        plot(x,e, 'r. ');%画出误差曲线
        xlabel('x轴'),ylabel('误差');
        title('步长为0.05时的误差曲线');

    大神,我修改完,还是运行不出来啊。你能不能直接把你的完整代码发给我,我自己找找错误。

  • hzlhm

    怎么发给你?

猜你喜欢