zoukankan      html  css  js  c++  java
  • 微分方程求解

    1. 欧拉方法

    通过计算:

      

    求[a,b]上的初值问题 y'=f(t,y),y(a)=y0的近似解。

    function E=euler(f,a,b,ya,M)
    %Input  - f is the function entered  as a string 'f'
    %       - a and b are the left and right end points
    %       - ya is the initial condition y(a)
    %       - M is the number of steps
    %Output - E=[T' Y'] where T is the vector of abscissas and
    %         Y is the vector of ordinates
    h=(b-a)/M;
    Y=zeros(1,M+1);
    T=a:h:b;
    Y(1)=ya;
    for j=1:M
        Y(j+1)=Y(j+h*feval(f,T(j),Y(j)));
    end
    E=[T' Y'];
    end

     2.四阶龙格-库塔法

    计算公式为:

      

    求[a,b]上的初值问题 y'=f(t,y),y(a)=y0的近似解。

    function R=rk4(f,a,b,ya,M)
    %Input  - f is the function entered  as a string 'f'
    %       - a and b are the left and right end points
    %       - ya is the initial condition y(a)
    %       - M is the number of steps
    %Output - R=[T' Y'] where T is the vector of abscissas and
    %         Y is the vector of ordinates
    %
    h=(b-a)/M;
    T=zeros(1,M+1);
    U=zeros(4,M+1);
    T=a:h:b; 
    Y(1)=ya;
    for j=1:M
        k1=h*feval(f,T(j),Y(j));
        k2=h*feval(f,T(j)+h/2,Y(j)+k1/2);
        k3=h*feval(f,T(j)+h/2,Y(j)+k2/2);
        k4=h*feval(f,T(j)+h,Y(j)+k3);
        Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)/6;
    end
    R=[T' Y'];
    end
  • 相关阅读:
    美化盒子和文本字体
    图片和多媒体
    学习node1_module对象
    学习vue5_组件
    学习vue4_input
    学习vue3
    学习vue2
    Ubuntu中U盘识别不了
    docker 建立新用户软件安装环境ubuntu
    计算机性能优化笔记
  • 原文地址:https://www.cnblogs.com/guojun-junguo/p/10139104.html
Copyright © 2011-2022 走看看