李拾梧 2022-06-01 10:50 采纳率: 100%
浏览 181
已结题

Matlab欧拉法求解非线性常微分方程

请教一下
这个怎么用Matlab编程求解这个非线性常微分方程
包括显式欧拉法,隐式欧拉法,两步欧拉法,改进欧拉法
包括图像

img

  • 写回答

1条回答 默认 最新

  • 一大岐 2022-06-01 13:20
    关注
    clear
    clc
    close all 
    %%  仿真步长 h=0.01 时
    Hfun = @(t,x) [ x(1)-x(2)-x(1)*x(3);x(1)+x(2)-x(2)*x(3); x(1)*x(1)+x(2)*x(2)-x(3)];   % 一阶微分方程导数表达式 
    % 参数
    t = [0 12];     % 时间范围
    h = 0.01;       % 时间步长
    x0 = [2 0 1];   % 初始状态 
    % 显式欧拉法求解
    [T1,X1] = ODE_ExplicitEuler( Hfun,t,h,x0 );
    % 隐式欧拉法求解
    [T2,X2] = ODE_ImplicitEuler( Hfun,t,h,x0 );
    % 两步欧拉法求解
    [T3,X3] = ODE_2OrderEuler( Hfun,t,h,x0 );
    % 改进欧拉法求解
    [T4,X4] = ODE_ImprovedEuler( Hfun,t,h,x0 );
    % Matlab自带ode45求解
    [T5,X5] = ode45( Hfun,t,x0 ); 
    % 绘图对比
    figure
    subplot(311)
    plot(T1,X1(:,1),T2,X2(:,1),T3,X3(:,1),T4,X4(:,1),T5,X5(:,1))
    xlabel('Time(s)')
    ylabel('x_1')
    legend('显式欧拉法','隐式欧拉法','两步欧拉法','改进欧拉法','ode45')
    subplot(312)
    plot(T1,X1(:,2),T2,X2(:,2),T3,X3(:,2),T4,X4(:,2),T5,X5(:,2))
    xlabel('Time(s)')
    ylabel('x_2')
    legend('显式欧拉法','隐式欧拉法','两步欧拉法','改进欧拉法','ode45')
    subplot(313)
    plot(T1,X1(:,3),T2,X2(:,3),T3,X3(:,3),T4,X4(:,3),T5,X5(:,3))
    xlabel('Time(s)')
    ylabel('x_3')
    legend('显式欧拉法','隐式欧拉法','两步欧拉法','改进欧拉法','ode45')
    % 显示欧拉
    function [T,X,dX] = ODE_ExplicitEuler( Hfun,t,h,x0 )
        % 确定时间节点
        T = t(1):h:t(2);
        % 计算
        N = length(T);
        x0 = x0(:);  x0 = x0';     % 初值变为行向量  
        m = length(x0);            % 状态量维数
        X = zeros(N,m);            % 初始化状态量
        dX = zeros(N,m);           % 状态导数
        X(1,:) = x0;
        for k = 2:N
            dX(k-1,:) = Hfun( T(k-1),X(k-1,:) );   
            h = T(k) - T(k-1);
            X(k,:) = X(k-1,:) + h*dX(k-1,:);                        
        end
        dX(N,:) = Hfun( T(N),X(N,:) );
    end
    % 隐式欧拉
    function [T,X,dX] = ODE_ImplicitEuler( Hfun,t,h,x0 )
        T = t(1):h:t(2);
        % 计算
        N = length(T);
        x0 = x0(:);  x0 = x0';     % 初值变为行向量  
        m = length(x0);            % 状态量维数
        X = zeros(N,m);            % 初始化状态量
        dX = zeros(N,m);           % 状态导数
        X(1,:) = x0;
        for k = 2:N
            h = T(k) - T(k-1);
            for jj = 1:m
               Ind =  zeros(1,m);         % 状态选择向量
               Ind(jj) = 1;        % 选择的向量位置  
               fh = @(x) h*Ind*Hfun(T(k),[X(k-1,1:jj-1)';x; X(k-1,jj+1:m)']) + X(k-1,jj)' - x;    
               X(k,jj) = fzero( fh,X(k-1,jj) ) ;   
            end                  
        end
        dX(N,:) = Hfun( T(N),X(N,:) );
    end
    
    function [T,X,dX] = ODE_2OrderEuler( Hfun,t,h,x0 )
        T = t(1):h:t(2);
        % 计算
        N = length(T);
        x0 = x0(:);  x0 = x0';     % 初值变为行向量  
        m = length(x0);            % 状态量维数
        X = zeros(N,m);            % 初始化状态量
        dX = zeros(N,m);           % 状态导数
        X(1,:) = x0;
        for k = 2:N
            dX(k-1,:) = Hfun( T(k-1),X(k-1,:) );   
            h = T(k) - T(k-1);
            if k == 2
                X(k,:) = X(k-1,:) + h*dX(k-1,:);    
            else
                X(k,:) = X(k-2,:) + 2*h*dX(k-1,:);   
            end
        end
        dX(N,:) = Hfun( T(N),X(N,:) );
    end
    function [T,X,dX] = ODE_ImprovedEuler( Hfun,t,h,x0 )
        T = t(1):h:t(2);
        % 计算
        N = length(T);
        x0 = x0(:);  x0 = x0';     % 初值变为行向量  
        m = length(x0);            % 状态量维数
        X = zeros(N,m);            % 初始化状态量
        dX = zeros(N,m);           % 状态导数
        X(1,:) = x0;
        for k = 2:N
            dX(k-1,:) = Hfun( T(k-1),X(k-1,:) );   
            h = T(k) - T(k-1);
            Xp = X(k-1,:) + h*dX(k-1,:);
            dXp = Hfun( T(k),Xp );
            X(k,:) = X(k-1,:) + (h/2)*(dX(k-1,:)+dXp');
        end
        dX(N,:) = Hfun( T(N),X(N,:) );
    end
    
    本回答被题主选为最佳回答 , 对您是否有帮助呢?
    评论 编辑记录

报告相同问题?

问题事件

  • 系统已结题 6月9日
  • 已采纳回答 6月1日
  • 创建了问题 6月1日

悬赏问题

  • ¥15 python使用pulp线性优化时报错
  • ¥15 开源或低价数据中台哪个最好
  • ¥15 arduino编程出现字符串疑似覆盖现象
  • ¥15 我的b站在没有碰到屏幕的情况下偶尔会自动跳出进度条,就像在屏幕上点了一下一样,但我并没有点。而且视频进度并没有变。这可能是什么原因造成的?
  • ¥30 STK matlab python仿真
  • ¥15 关于IMageEnView 图标定位问题
  • ¥20 求解答(matlab)
  • ¥30 ffmpeg库使用过程中遇到的问题
  • ¥15 pyqt5 中python如何通过Qtwebchannel主动发消息给web前端
  • ¥15 关于HTML中title获取xml内容的问题