Ode23，ode45和ode113在Scilab

Question

我在Scilab中使用ode23、ode45和ode113来求解ODE，我将计算平均绝对误差来检验哪一个更准确，但是所有的求解者都给出误差等于0，向量长度对所有的n=100都是相同的。我看到这里有一些问题，它不符合MATLAB的结果。你能帮忙吗？我要附上我的密码。

clear
function [v] = functiondifference(t,y,exact)

    [row,column] = size(t); // size(A) returns the dimensions of a matrix A. In our case t is a column vector
    for k = 1:row
       v(k) = feval(t(k),exact) - y(k);   /////////////////////note
    end
endfunction
function y = exact(t)

  y = -3-exp(-sin(t)/2);

endfunction
function yp=de(t,y)
    yp=-(3+y/2)*cos(t)
endfunction

function dydt=f(t,y)
    y=y(1);
    dydt=[ -(3+y/2)*cos(t)]
endfunction
t=linspace(0,%pi);
y0=-4;
//ode 23
y = ode("adams", y0, 0, t, f);   //t0=0
err=functiondifference(t,y,exact)
Error=mean(abs(err))
L=length(y)
disp(L)
//ode 45
y1 = ode("rkf", y0, 0, t, f);
err1=functiondifference(t,y,exact)
Error1=mean(abs(err1))
L1=length(y1)
disp(L1)
//ode 113
y2= ode(y0, 0, t, f);  
err2=functiondifference(t,y2,exact)
Error2=mean(abs(err2))
L2=length(y2)
disp(L2)

谢谢你的帮助。对于给定的微分方程yp=0.5*(3+y)*cos(t)，精确解是正确的，即使函数改为行，列=大小(t‘)，输出也是不正确的。我要附上

MATLAB的结果和精确解。提前谢谢。

Answer 1

错误/误解：

• ，您在注释中写到t是列向量。但实际上它是一个行向量。因此，您做了一个不同的计算，在初始值上没有差别。最简单的修正是，使用t的转置来使假设成为现实：

行，列=大小(t‘)；

• 您的确切解决方案不是ODE yp = -(3+y/2)*cos(t) = -0.5*(6+y)*cos(t)的确切解决方案。分离后，这会给

(6+y(0))*exp(-0.5*sin(t)) = 6+y(t)

通过这两种修正，我得到了adams、rkf和lsoda的错误1.953D-07, 4.954D-12, 5.978D-07。

根据注释，为了获得更好的洞察力，增强脚本来计数函数计算和测试多个错误公差。

clear
function y = exact(t)
  y = -6+2*exp(-sin(t)/2);
end//function

global fcnt;

function yp=de(t,y)
    global fcnt
    fcnt = fcnt + 1;
    yp = -(3+y/2)*cos(t)
end//function

t=linspace(0,%pi,10);
y0=exact(t(1));

method=["lsode","adams","rkf","stiff"]
for k = 1:4
    for j=5:8
        tol = 10^(-j);
        fcnt = 0;
        if k==1
            y = ode(y0, 0, t, 0.1*tol, tol, de); 
        else
            y = ode(method(k), y0, 0, t, 0.1*tol, tol, de);   //t0=0
        end
        err=y-exact(t);
        Error=mean(abs(err));//*t(length(t));
        printf("%6s: toll=%6.3g, err=%6.3g, fcnt=%6d\n",method(k),tol,Error,fcnt)
    end
end

这在表中得到了结果。

 lsode: toll= 1e-05, err=2.32e-05, fcnt=    57
 lsode: toll= 1e-06, err=5.81e-07, fcnt=    63
 lsode: toll= 1e-07, err=8.43e-08, fcnt=    81
 lsode: toll= 1e-08, err=6.41e-08, fcnt=   107
 adams: toll= 1e-05, err=2.21e-05, fcnt=    40
 adams: toll= 1e-06, err=5.88e-07, fcnt=    47
 adams: toll= 1e-07, err=7.45e-08, fcnt=    58
 adams: toll= 1e-08, err=4.29e-08, fcnt=    72
   rkf: toll= 1e-05, err=7.35e-07, fcnt=    61
   rkf: toll= 1e-06, err=7.39e-07, fcnt=    61
   rkf: toll= 1e-07, err=2.24e-07, fcnt=   113
   rkf: toll= 1e-08, err=2.35e-08, fcnt=   132
 stiff: toll= 1e-05, err=2.70e-05, fcnt=    72
 stiff: toll= 1e-06, err=3.40e-06, fcnt=    60
 stiff: toll= 1e-07, err=5.09e-07, fcnt=    87
 stiff: toll= 1e-08, err=1.55e-07, fcnt=   121

• 方法满足精度目标
• “RKF”不采用阶内插补，因此隐式地将输出样本的步长作为内步长的上限。这导致比necessary.
• "adams“更多的函数调用确实是最便宜的方法，但是作为显式求解器
• ，情景隐式求解器”默认/lsode“和”刚性“在中间有一个相似的轮廓，显然隐式步骤允许稍微大一点的步长，否则就会期望显式多步方法”adams“将计数加倍。

Answer 2

给定DE的精确解

​