Matlab中的ode45中变参数迭代

写在前面

ode45是Matlab中处理微分方程中常用的求解函数，但是在实际解决问题中，方程总是多种多样的，并不是普通形式的。比如我遇到的问题，右边的微分方程中的系数是变化的，不是常数，同时并不是根据横坐标x变化，而是根据函数值y进行变化的。

ode45积分器

这里我们先求解一个最简单的微分方程：

y'=2t

积分区间[0,5],初始条件y0 = 0。

tspan = [0 5];
y0    = 0;
[t,y] = ode45(@(t,y) 2*t, tspan,y0);
%%画图
plot(t,y,'-o')

结果： 

解算非刚性方程

我们有如下二阶方程

这里做代换，用y2作为y1的导数，y1‘=y2，可重写为一阶ODE方程组。

通过调用函数可以实现多个导数的求解，代码如下：

[t,y] = ode45(@func1,[0 20],[2,0]);
plot(t,y(:,1),'-o',t,y(:,2),'-o')
title('Solution of van der Pol Equation (\mu = 1) with ODE45');
xlabel('Time t');
ylabel('Solution y');
legend('y_1','y_2')

function dydt = func1(t,y)
dydt = zeros(2,1);
dydt(1) = y(2);
dydt(2) = (1-y(1)^2)*y(2)-y(1);
end

上述案例我们通过调用函数解决了多个（高阶）微分未知量。

随y变化的微分方程

现在我们有方程

积分区间为[1,5],初始值y(0) = 1.

代码：

tspan = [1 5]
y0    = 1
[t,y] = ode45(@(t,y) func1(t,y,i),tspan,y0)
plot(t,y,'-o')

function dydt = func1(t,y,i)
k    = y+1;
dydt = k+2*y+1;
end

结果

但是这里有个问题是，每次迭代使用的y值其实不是最后输出的y值，两者之间存在误差。

我们更改代码查看ode45内部的参数变化，代码修改为：

tspan = [1 5]
y0    = 1
[t,y] = ode45(@(t,y) func1(t,y,i),tspan,y0)
plot(t,y,'-o')

function dydt = func1(t,y,i)
k = y+1；
y
dydt = k+2*y+1;
end

左侧为函数运行中的y值，右侧为输出的y值。

如左侧的第5个数为1.1892，右侧的第5个数为1.2136。两者并不相同，作者也不明白两者存在误差的原因。

如果有大佬知道其中计算的过程以及误差的问题，欢迎在评论区或者私信笔者进行讨论。