![]() |
哪位大侠能帮小弟改一下这个有关三次样条 插值多项式的程序
w=imread('11.jpg');
[x,y]=size(w); n = size(x,2); % 插值点个数, 注意: 不是书上的 n for j = 1 : n-1 h(j) = x(j+1) - x(j); % 步长 f(j) = (y(j+1) - y(j)) / h(j); % 一阶差商 end for j = 2 : n-1 lambda(j) = h(j)/(h(j-1) + h(j)); mu(j) = h(j-1)/(h(j-1) + h(j)); g(j) = 3 * (lambda(j) * f(j-1) + mu(j) * f(j)); end A = 2 * eye(n,n); for i = 2 : n-1 A(i,i+1) = mu(i); A(i,i-1) = lambda(i); end b = zeros(n,1); for i = 2 : n-1 b(i) = g(i); end % 2 阶导数边界条件 %b(1) = 3*f(1); %A(1,2) = 1; %b(n) = 3*f(n-1); %A(n,n-1) = 1; m = A \ b; % 求解线性方程组, 得到各节点的一阶导数值 S = []; for j = 1 : n-1 s = (y(j)/(h(j)^3))*conv([1 -2*x(j+1) x(j+1)^2],[2 h(j)-2*x(j)]) + ... (y(j+1)/(h(j)^3))*conv([1 -2*x(j) x(j)^2],[-2 h(j)+2*x(j+1)]) + ... (m(j)/(h(j)^2))*conv([1 -2*x(j+1) x(j+1)^2],[1 -x(j)]) + ... (m(j+1)/(h(j)^2))*conv([1 -2*x(j) x(j)^2],[1 -x(j+1)]); S = [S;s]; end |
所有时间均为北京时间。现在的时间是 05:26。 |
Powered by vBulletin
版权所有 ©2000 - 2025,Jelsoft Enterprises Ltd.