中点法求解高阶线性微分方程（Solving higher order linear differential equations by midpoint method）

欧拉法（Euler）和改进欧拉法是求解线性微分方程初值问题常用的数值解法，对于一般方程可以手算出来，但是对于比较复杂的手算是算不出来，本文以二阶方程为例编写matlab程序实现高阶微分方程的数值解法。并且估计了绝对误差与时间步长h之间的关系。

%% 定义计算步长，设置变量的初始值
clear;clc;close all
Delta=0.001;% 定义步长
t=0:Delta:1;% 定义自变量t
n=length(t);
Y(:,1)=[1;0] % 定义x y的初始值
%% 自定义改进欧拉法，求解微分方程组
for k=1:n-1
    F1=f(t(k),Y(:,k));
    F2=f(t(k+1),Y(:,k)+Delta*f(t(k),Y(:,k)));
    Y(:,k+1)=Y(:,k)+Delta*(F1+F2)/2;
end
x=Y(1,:);
y=Y(2,:);
%% 绘制x y的求解结果
figure
set(gcf,'units','normalized','position',[0.15 0.2 0.7 0.6]); %设置figure窗口的位置和尺寸
subplot(1,2,1)
plot(t,x,'b')
xlabel('t')
ylabel('x')

subplot(1,2,2)
plot(t,y,'r')
xlabel('t')
ylabel('y')

求解误差随时间步长的变化

%% 求误差值
clear;
clc;
err1=0;
err2=0;
j=1;
for h=0.001:0.001:0.1
    N=1/h;
    y(:,1)=[1
            0];
    a=[0  , 1
       -1 , -1];
  for i = 1:N
    t=i*h;
    b1=[0
        sin(t)];
    b2=[0
        sin(t+h/2)];
    f1=y(:,i)+1/2*h*(a*y(:,i)+b1); 
    y(:,i+1) = y(:,i) + h*(a*f1+b2);
    err1=err1+h^3/12*abs(y(1,i)+cos(t)-sin(t));
    err2=err2+h^3/12*abs(y(2,i)-sin(t)-cos(t));
  end
    err(j)=sqrt(err1^2+err2^2);
    h1(j)=h;
    j=j+1;
end
 %% 绘制误差图
figure
set(gcf,'units','normalized','position',[0.15 0.2 0.7 0.6]); %设置figure窗口的位置和尺寸
subplot(1,2,1)
plot(h1,err,'b')
xlabel('h1')
ylabel('err')

subplot(1,2,2)
plot(log(h1),log(err),'r')
xlabel('log(h1)')
ylabel('log(err)')