Matlab - 微分方程求解器 Ode45 的意外结果

Question

我正在尝试使用 MATLAB 的 ode 求解器 ode45 求解一个微分方程。我已经尝试将它与其他更简单的函数一起使用，并让它绘制函数。它们看起来都正确，但是当我插入需要解决的功能时，它会失败。该图从 y(0) = 1 开始，但在某个点开始下降，此时它应该一直是一个递增函数，一直到其临界点。

function [xpts,soln] = diffsolver(p1x,p2x,p3x,p1rr,y0)
syms x y
yp = matlabFunction((p3x/p1x) - (p2x/p1x) * y); 
[xpts,soln] = ode45(yp,[0 p1rr],y0); 

p1x、p2x 和 p3x 是多项式，它们作为参数传递到此 diffsolver 函数中。

p1rr 这里是关键点。函数应该在临界点之后发散，所以我想把它整合到那个点。

编辑：这是我在使用 diffsolver 之前的代码，即上述函数。我做了填充近似来找到多项式 p1、p2 和 p3。然后我找到临界点，即最接近目标的 p1 的根（目标由用户指定）。

我检查临界点是否为空（有时某些函数中可能没有临界点）。如果不为空，则使用上述函数求解微分方程。然后它基本上绘制了从上述函数返回的 x 和 y 点。

function error = padeapprox(m,n,j)
global f df p1 p2 p3 N target

error = 0;
size = m + n + j + 2;

A = zeros(size,size);
for i = 1:m
    A((i + 1):size,i) = df(1:(size - i));
end
for i = (m + 1):(m + n + 1)
    A((i - m):size,i) = f(1:(size + 1 - i + m));
end
for i = (m + n + 2):size
    A(i - (m + n + 1),i) = -1;
end

if det(A) == 0
    error = 1;
    fprintf('Warning: Matrix is singular.\n');
end

V = -A\df(1:size); 

p1 = [1];
for i = 1:m
    p1 = [p1; V(i)];
end

p2 = [];
for i = (m + 1):(m + n + 1)
    p2 = [p2; V(i)];
end

p3 = [];
for i = (m + n + 2):size
    p3 = [p3; V(i)];
end

fx = poly2sym(f(end:-1:1)); 
dfx = poly2sym(df(end:-1:1));
p1x = poly2sym(p1(end:-1:1)); 
p2x = poly2sym(p2(end:-1:1)); 
p3x = poly2sym(p3(end:-1:1)); 
p3fullx = p1x * dfx + p2x * fx; 
p3full = sym2poly(p3fullx); p3full = p3full(end:-1:1); 

p1r = roots(p1(end:-1:1));
p1rr = findroots(p1r,target); % findroots eliminates unreal roots and chooses the one closest to the target

if ~isempty(p1rr)
    [xpts,soln] = diffsolver(p1x,p2x,p3fullx,p1rr,f(1));
    if rcond(A) >= 1e-10
        plot(xpts,soln); axis([0 p1rr 0 5]); hold all
    end
end

我看到了一些使用另一个函数来生成微分方程的示例，但我尝试过将 matlabFunction() 方法与其他更简单的函数一起使用，它似乎可以工作。只是当我尝试解决此功能时，它失败了。当解决的值都应该为正时，它们开始变为负。

我还尝试使用另一个求解器 dsolve()。但它一直给我一个隐含的解决方案......

有人知道为什么会这样吗？任何建议表示赞赏。谢谢！

Answer 1

由于您的代码似乎适用于更简单的功能，您可以尝试增加 ode45 求解器的精度选项。

这可以通过使用 odeset 来实现：

 options = odeset('RelTol',1e-10,'AbsTol',1e-10);
 [T,Y] = ode45(@function,[tspan],[y0],options);