写了一个自适应 Simpson 求积公式的代码, 运行过程中出现了 abs命令或者变量的类型发生变化导致,无法求值 的错误, 如图所示我写的 adapsimp 的代码如下
function [s, err] = adapsimp(func, a, b, tol) s = comsimp(func, a, b, 2); c = (a + b) / 2; s1 = comsimp(func, a, c, 2); s2 = comsimp(func, c, b, 2); s12 = s1 + s2; err = abs(s12 - s) / 15; if err < tol s = s12; else [s1, err1] = adapsimp(func, a, c, tol/2); [s2, err2] = adapsimp(func, c, b, tol/2); s = s1 + s2; err = err1 + err2; end end里面用到了 comsimp 函数, 是这样写的
function s = comsimp(func, a, b, n) h = (b - a) / n; s0 = func(a) + func(b); s1 = 0; % summation of f(x_{2k-1}) s2 = 0; % summation of f(x_{2k}) for k = 1:n-1 x = a + k * h; if rem(k , 2) == 0 s2 = s2 + func(x); else s1 = s1 + func(x); end end s = h * (s0 + 4 * s1 + 2 * s2) / 3; end这部分代码在 octave 上运行是没有问题的用的版本是2.2.0最新版的.
h=0.01; r=0.5; tao=r*h; x=0:h:5; t=0:tao:10; length_t=length(t); length_x=length(x); U=zeros(length_t,length_x); length_in=length_x-1; A=diag((1-r)*ones(length_in,1))+r*diag(ones(length_in-1,1),-1); [X,T]=meshgrid(x,t); F=2*T.*sin(X)+T.^2.*cos(X); for i=2:length(t) U(i,2:end)=A*U(i-1,2:end)'+tao*F(i,2:end)'; end U_acc=T.^2.*sin(X); figure(1) subplot(1,2,1) mesh(x,t,U) title("解曲面") subplot(1,2,2) mesh(x,t,U_acc-U) title("误差曲面")发现是FOR循环内的矩阵转置处出现了问题,但是无法清晰定位具体发生了什么导致的问题,希望开发组能解决。