Matlab各种插值法.part6.zip_newton_三次样条插值_二维插值_插值_插值matlab

function [ybar,ybarSym ] = Fredholm2(ker,f,lambda,a,b,n) %求解第二类Fredholm积分方程:y(x) = lambda*int(k(x,s)*y(s)*ds,a,b)+f(x),a<=x<=b %ker,核函数k(x,s)的函数句柄；f，自由项f(x)的函数句柄；lambda，参数；a,b积分上下限 %n,Gauss-Lobatto积分公式的积分节点数 % ybar，y(x)的近似函数，以函数句柄形式返回 % ybarSym，y(x)的近似函数，与ybar等价，为方便查看以符号公式形式返回 eval(['ker=',vectorize(ker),';']);%ker统一为点乘形式 eval(['f=',vectorize(f),';']);%f统一为点乘形式 if nargin<=5 %没有指定n的话，n为5 n = 5; end P = LegendrePoly(n);%调用LegendrePoly求n次Legendre多项式 diffP = diff(P);%n次n次Legendre多项式的导数 %上面得到的diffP还是符号多项式表达式，利用sym2poly函数将其转化成一元n次多项式的标准 %形式，即从最高次到常数项的所有系数组成的向量，利用roots求其根，并利用sort从小到大排 %好序和-1，1一起得到n次Legendre多项式的求积节点 xk = [-1;sort(roots(sym2poly(diffP)));1]; %计算求积系数 Ak = (2/(n^2+n)).*[1;1./(polyval(sym2poly(P),xk(2:end-1)).^2);1]; if a~=-1 || b~=1 %如果[a,b]不是[-1,1]的话，进行变量替换 tmp = (b-a)/2; xab = tmp*xk + (a+b)/2; F = f(xab);%得到自由项f(x)关于求积节点的函数值向量 if isscalar(F)%针对某些自由项函数为常数的情形 F = F*ones(size(xab)); end [Y1,X1] = meshgrid(xab,xab); kerXY = ker(X1,Y1); if isscalar(kerXY)%针对某些核函数为常数的情形 kerXY = kerXY*ones(size(X1)); end K = bsxfun(@times,kerXY,Ak');%y的线性方程组的系数矩阵 K = (b-a)/2*lambda*K;%进行积分变量替换后的K Y = (eye(n+1)-K)\F;%求解线性方程组得到Y X = xab;%对应的自变量X值 else %积分区间为[-1,1]，不用变量替换 F = f(xk); if isscalar(F) F = F*ones(size(xk)); end [Y1,X1] = meshgrid(xk,xk); kerXY = ker(X1,Y1); if isscalar(kerXY)%针对某些核函数为常数的情形 kerXY = kerXY*ones(size(X1)); end K = lambda*bsxfun(@times,kerXY,Ak');%y的线性方程组的系数矩阵 Y = (eye(n+1)-K)\F; X = xk; end fun = cell(n+1,1); for k=1:n+1 %积分离散化后每一个节点上的分量函数 fun{k} = @(x) Ak(k)*Y(k)*ker(x,X(k)); end ybar = @approxY;%返回y(x)的近似函数 function yx = approxY(x) %利用嵌套函数形式可以方便的计算ybar的近似函数 subY = fun{1}(x); for jj=2:n+1 subY = subY + fun{jj}(x); end yx = f(x)+ (b-a)/2*lambda*subY; end %下面利用符号计算得到ybar的数学表达式，读者可以根据需要采用该方法或者采用 %上面利用嵌套函数返回函数句柄的做法 ybarSym = GenerateExpr(f,a,b,lambda,Ak,ker,X,Y); end %======================================== %计算n次Legendre多项式的函数(普通子函数形式) %======================================== function P = LegendrePoly(n) syms x; if n == 1 P = sym(1); return else P = simplify( diff((x^2-1)^n,n)/2^n/gamma(n+1)); end end %======================================================= %符号计算得到ybar的数学表达式,之所以采用普通子函数是因为 %含有嵌套函数情况下，主函数或者嵌套函数不允许动态生成符号变量 %======================================================== function ybarSym = GenerateExpr(f,a,b,lambda,Ak,ker,X,Y) digits(16); As = sym(Ak,'d'); Xs = sym(X,'d'); Ys = sym(Y,'d'); syms x; ybarSym = f(x)+ (b-a)/2*lambda*sum(As.*ker(x,Xs).*Ys); ybarSym = simplify(ybarSym); end