NEWTON插值处理，包含主函数、两个子函数

% 待求解方程 function y = fun(x) format long; global ca cb ka kb kw n = size(kb,2); m = size(ka,2); %% 方程左半部分 % 方程左侧第二项分母项 if n == 1 term12 = (kw/x) + kb(1); elseif n==2 term12 = (kw/x)^2 + kb(1)*(kw/x) + kb(1)*kb(2); elseif n==3 term12 = (kw/x)^3 + kb(1)*(kw/x)^2 + kb(1)*kb(2)*(kw/x) + kb(1)*kb(2)*kb(3); elseif n==4 term12 = (kw/x)^4 + kb(1)*(kw/x)^3 + kb(1)*kb(2)*(kw/x)^2 + kb(1)*kb(2)*kb(3)*(kw/x) + kb(1)*kb(2)*kb(3)*kb(4); elseif n==5 term12 = (kw/x)^5 + kb(1)*(kw/x)^4 + kb(1)*kb(2)*(kw/x)^3 + kb(1)*kb(2)*kb(3)*(kw/x)^2 + kb(1)*kb(2)*kb(3)*kb(4)*(kw/x) + kb(1)*kb(2)*kb(3)*kb(4)*kb(5); elseif n==6 term12 = (kw/x)^6 + kb(1)*(kw/x)^5 + kb(1)*kb(2)*(kw/x)^4 + kb(1)*kb(2)*kb(3)*(kw/x)^3 + kb(1)*kb(2)*kb(3)*kb(4)*(kw/x)^2 + kb(1)*kb(2)*kb(3)*kb(4)*kb(5)*(kw/x) + kb(1)*kb(2)*kb(3)*kb(4)*kb(5)*kb(6); end term1 = 0; for i = 1:n if i == 1 term11 = i*kb(1)*(kw/x)^(n-i); % 方程左侧第二项分子项 elseif i==2 term11 = i*kb(1)*kb(2)*(kw/x)^(n-i); elseif i==3 term11 = i*kb(1)*kb(2)*kb(3)*(kw/x)^(n-i); elseif i==4 term11 = i*kb(1)*kb(2)*kb(3)*kb(4)*(kw/x)^(n-i); elseif i==5 term11 = i*kb(1)*kb(2)*kb(3)*kb(4)*kb(5)*(kw/x)^(n-i); elseif i==6 term11 = i*kb(1)*kb(2)*kb(3)*kb(4)*kb(5)*kb(6)*(kw/x)^(n-i); end term1 = (term1 + term11/term12); %方程左侧第二项 end %% 方程右半部分 % 方程右侧第二项分母项 if m == 1 term22 = x + ka(1); elseif m==2 term22 = (x)^2 + ka(1)*(x) + ka(1)*ka(2); elseif m==3 term22 = (x)^3 + ka(1)*(x)^2 + ka(1)*ka(2)*(x) + ka(1)*ka(2)*ka(3); elseif m==4 term22 = (x)^4 + ka(1)*(x)^3 + ka(1)*ka(2)*(x)^2 + ka(1)*ka(2)*ka(3)*(x) + ka(1)*ka(2)*ka(3)*ka(4); elseif m==5 term22 = (x)^5 + ka(1)*(x)^4 + ka(1)*ka(2)*(x)^3 + ka(1)*ka(2)*ka(3)*(x)^2 + ka(1)*ka(2)*ka(3)*ka(4)*(x) + ka(1)*ka(2)*ka(3)*ka(4)*ka(5); elseif m==6 term22 = (x)^6 + ka(1)*(x)^5 + ka(1)*ka(2)*(x)^4 + ka(1)*ka(2)*ka(3)*(x)^3 + ka(1)*ka(2)*ka(3)*ka(4)*(x)^2 + ka(1)*ka(2)*ka(3)*ka(4)*ka(5)*(x) + ka(1)*ka(2)*ka(3)*ka(4)*ka(5)*ka(6); end term2 = 0; for i = 1:m if i == 1 term21 = i*ka(1)*x^(m-i); % 方程右侧第二项分子项 elseif i==2 term21 = i*ka(1)*ka(2)*x^(m-i); elseif i==3 term21 = i*ka(1)*ka(2)*ka(3)*x^(m-i); elseif i==4 term21 = i*ka(1)*ka(2)*ka(3)*ka(4)*x^(m-i); elseif i==5 term21 = i*ka(1)*ka(2)*ka(3)*ka(4)*ka(5)*x^(m-i); elseif i==6 term21 = i*ka(1)*ka(2)*ka(3)*ka(4)*ka(5)*ka(6)*x^(m-i); end term2 = (term2 + term21/term22); % 方程右侧第二项 end y = x + cb*term1 - kw/x - ca*term2; %待求方程