微分方程的数值解法与符号解法

数值解

问题

Matlab实现

main.m 主函数

clear,clc
x0 = 0;  % xmin
x1 = 1;  % xmax
h = 0.1;  % 步长
y0 = 1/2;   % y0为初始值

% 欧拉法
[xy{1}(:,1), xy{1}(:,2)] = forwardEuler(x0, y0, x1, h);

% 改进欧拉法
[xy{2}(:,1), xy{2}(:,2)] = improvedEuler(x0, x1, y0, h);

% 四阶龙格库塔法
[xy{3}(:,1), xy{3}(:,2)] = runge(x0, x1, y0, h);

% 解析解
% p(x) = -3;
% q(x) = x;  
f = @(x) (11*exp(3*x))/18 - x/3 - 1/9;
x1 = 0:0.1:1;
y1 = f(x1);

% 画图
hold on       
plot(xy{1}(:,1), xy{1}(:,2),'-o','LineWidth', 1.5)
plot(xy{2}(:,1), xy{2}(:,2),'-d','LineWidth', 1.5)
plot(xy{3}(:,1), xy{3}(:,2),'-s','LineWidth', 1.5)
plot(x1,y1,'linewidth',1.5)

legend('欧拉法','改进欧拉法','四阶龙格库塔法','解析解','location','best')

% 计算误差
for i = 1:3
    error(:,i) = y1' - xy{i}(:,2);   % 截断误差
end
error2 = abs(error);   % 误差绝对值

func.m 微分方程函数

%func.m文件  常微分方程为y' = x+3y
function z = func(x, y)
z = x + 3*y;

forwardEuler.m前向欧拉法/欧拉法

%forwardEuler.m文件
function [x, y] = forwardEuler(x0, y0, x1, h)
n = floor((x1 - x0)/h);
x = zeros(n + 1, 1);
y = zeros(n + 1, 1);
x(1) = x0;
y(1) = y0;
for i = 1 : n
    x(i + 1) = x(i) + h;
    y(i + 1) = y(i) + h * func(x(i), y(i));
end

improvedEuler.m改进欧拉法

%改进的欧拉算法代码
function [x, y] = improvedEuler(x0, x1, y0, h)
n = floor((x1 - x0)/h);
x = zeros(n + 1, 1);
y = zeros(n + 1, 1);
x(1) = x0;
y(1) = y0;
for i = 2 : n+1
    x(i) = x(i-1) + h;
    yp = y(i-1) + h*func(x(i-1),y(i-1));
    yc = y(i-1) + h*func(x(i),yp);
    y(i) = (yp+yc)/2;
%     y(i) = y(i-1) + h * (y(i - 1) - 2 * x(i - 1) / y(i - 1));
%     % 用上一步求得的值代入这一步，变为显式算法
%     y(i) = y(i - 1) + h / 2 * (y(i - 1) - 2 * x(i - 1) / y(i - 1) + y(i) - 2 * x(i) / y(i));
end

runge.m龙格库塔法

function [x, y] = runge(x0, x1, y0, h)
n = (x1 - x0) / h;
x = zeros(n + 1,1);
y = zeros(n + 1,1);
x(1) = x0;
y(1) = y0;
for i = 1:n
    x(i + 1) = x(i) + h;
    k1 = func(x(i), y(i));
    k2 = func(x(i) + 0.5*h, y(i) + k1*h/2);
    k3 = func(x(i) + 0.5*h, y(i) + k2*h/2);
    k4 = func(x(i)+ h, y(i) + k3*h);
    y(i + 1) = y(i) + h*(k1 + 2*k2 + 2*k3 + k4)/6;
end
end

结果

符号解

syms x y
dsolve('Dy=t+3*y','y(0) = 1/2')   % 解符号解