数值分析算法 MATLAB 实践 常微分方程求解

Euler 法及改进算法

function [x,y] = euler(fun,a,b,h,y0)
%一阶常微分方程的一般表达式的右端函数：fun
% 显示欧拉格式
% f是带求函数的一阶导形式
% a,b分别是自变量取值上下限
% y0 是初始条件y(0)
% h是步长
    s = (b - a) / h; % 求步数
    X = zeros(1, s+1);
    Y = zeros(1, s+1);
    X = a:h:b;
    Y(1) = y0;
    for k = 1:s
        Y(k+1) = Y(k) + h * fun(X(k), Y(k))
    end
    x = X';
    y = Y';
end

%% euler求解微分方程
% dfun1 = y^2-y^3;
 [x, y] = euler(@dfun1, 0,5,0.01,0.1);
% dfun1 = y;
% [x, y] = euler(@dfun1, 0,1,0.1,1);
figure
plot(x, y);
title('显示欧拉格式');

%% 微分方程
function dfun1 = dfun1(t,y)
    dfun1 = y^2-y^3;
    %dfun1 = y;
end

function[x,y]=imp_euler(func,a_start,b_end,h_step,y0)
%一阶常微分方程的一般表达式的右端函数：fun
% 显示欧拉格式
% func是带求函数的一阶导形式
% a_start,b_end分别是自变量取值上下限
% y0是初始条件y(0)
% h_step是步长
x = a_start : h_step : b_end;
N = length(x);
y = zeros(1, N);
y(1) = y0;
for i = 2:N
    % 显式 Euler 作为初始值迭代计算
    yi_0 =  y(i-1) + h_step * func(x(i-1), y(i-1));
    yi_1 = y(i - 1) + h_step * func(x(i), yi_0);
    while abs(yi_1 - yi_0) > 1e-6
        yi_0 = yi_1;
        yi_1 = y(i - 1) + h_step * func(x(i), yi_0);
    end
    y(i) = yi_1;
end

function[x,y]=improve_euler(func,a_start,b_end,h_step,y0)
%一阶常微分方程的一般表达式的右端函数：fun
% 显示欧拉格式
% func是带求函数的一阶导形式
% a_start,b_end分别是自变量取值上下限
% y0是初始条件y(0)
% h_step是步长
x = a_start : h_step : b_end;
N = length(x);
y = zeros(1, N);
y(1) = y0;
for i = 2:N
    yp = y(i-1) + h_step * func(x(i-1), y(i-1));
    yq = y(i-1) + h_step * func(x(i), yp);
    y(i) = 0.5 * (yp + yq);
end

Runge-Kutta 算法

4阶-单变量龙格库塔公式

4阶-多变量龙格库塔公式

% 单变量龙格库塔Runge_kutta 经典法
%一阶常微分方程的一般表达式的右端函数：func
% func是带求函数的一阶导形式
% a_start,b_end分别是自变量取值上下限
% y0是初始条件y(0)
% h_step是步长
function[x,y]=Runge_kutta(func,a_start,b_end,h_step,y0)
x = a_start : h_step : b_end;
N = length(x);
y = zeros(1, N);
y(1) = y0;
for i = 2:N
    k1 = func(x(i-1), y(i-1));
    k2 = func(x(i-1) + h_step/2, y(i-1) + h_step/2*k1);
    k3 = func(x(i-1) + h_step/2, y(i-1) + h_step/2*k2);
    k4 = func(x(i-1) + h_step, y(i-1) + h_step*k3);
    y(i) = y(i-1) + h_step/6*(k1 + 2*k2 + 2*k3 + k4);
end

function[x,y]=Runge_kutta45(dyfunc,xspan,y0,h_step)
% 多变量龙格库塔Runge_kutta45
% h_step是步长常选取为0.01；
% ufunc是函数名；
% x0是初始时间值；
% y0是初始化值； 
% n 是迭代步数；
      x = xspan(1):h_step:xspan(2);
      y = zeros(length(y0),length(x));
      y(:,1) = y0(:);
      %循环迭代数值求解部分
     for n=1 : (length(x)-1)
          k1=feval(dyfunc, x(n),y(:,n));
          k2=feval(dyfunc, x(n)+h_step/2,y(:,n)+h_step/2*k1);
          k3=feval(dyfunc, x(n)+h_step/2,y(:,n)+h_step/2*k2);
          k4=feval(dyfunc, x(n+1),y(:,n)+h_step*k3);
          y(:,n+1)=y(:,n)+h_step*(k1+2*k2+2*k3+k4)/6; 
          %按照4阶多变量龙格库塔方法进行数值求解
     end
end

clc;
clear all;
y0=[0,2,9];%初值
xspan = [0,200];%求解区间
h_step = 0.001;%ode45是变步长的算法
[x,y] = Runge_kutta45(@lorenz_diff,xspan,y0,h_step);
figure(1);
plot3(y(1,:),y(2,:),y(3,:),'.');title("x-y-z");
figure(2);
plot3(y(1,:),y(3,:),y(2,:),'.');title("x-z-y");
figure(3);
plot3(y(2,:),y(1,:),y(3,:),'.');title("y-x-z");
function dydt = lorenz_diff(t,y)
%{
    x-->y(1),y-->y(2),z-->y(3)
%}
dydt = [        10*(y(2)-y(1));
             -y(1)*y(3)+30*y(1)-y(2)
              y(1)*y(2)-8/3*y(3)];  
end

Matlab 函数库求解

[t, Xt] = ode45(odefun, tspan, X0)
odefun是函数句柄，可以是函数文件名，匿名函数句柄或内联函数名
tspan是区间 [t0 tfinal] 或者一系列散点[t0,t1,…,tf]
X0是初始值向量
t返回列向量的时间点
Xt返回对应T的求解列向量

Lorenz系统

function dydt = lorenz_diff(t,y)
%{
    x-->y(1),y-->y(2),z-->y(3)
%}
dydt = [        10*(y(2)-y(1));
             -y(1)*y(3)+30*y(1)-y(2)
              y(1)*y(2)-8/3*y(3)];  
end
clc;
y0 = [0,2,9];
[t,y] = ode45('lorenz_diff',[0,200],y0); 
%% 调用ode45绘制Lorenz系统 2D
figure(1);
plot(y(:,1),y(:,3),'.');
xlabel('x');ylabel('z');title("x-z");
figure(2);
plot(y(:,1),y(:,2),'.');
xlabel('x');ylabel('y');title("x-y");
figure(3);
plot(y(:,2),y(:,3),'.');
ylabel('y');zlabel('z');title("y-z");

%% 火焰传播数学模型求解
clc;
clear all;
delta=0.01;
f=@(t,y)y^2-y^3;
opts=odeset('Reltol',1.e-4);
[t1,y1]=ode45(f,[0  2/delta], delta, opts);
figure(1)
plot(t1,y1,'-','Marker','.');
title('数值解曲线');
ylabel('y'); xlabel('t');