《计算方法上机作业求三次样条插值函数的matlab程序》由会员分享，可在线阅读，更多相关《计算方法上机作业求三次样条插值函数的matlab程序（4页珍藏版）》请在金锄头文库上搜索。
1、计算方法上机报告 23 附录 3 求三次样条插值函数的 matlab 程序 clc;clear; type = input('n 选择插值节点的类型：n 插值对象为连续函数请输入 1，n 插值对象为列表函数请输入 2。n'); while type =1 pause(1) type = input('n 选择插值节点的类型：n 插值对象为连续函数请输入 1， n 插值对象为列表函数请输入 2。n'); end if type=1 s = input('n 请输入连续函数的表达式：nf(x) = ','s'); a =
2、input('n 请输入插值区间下限 a：n'); b = input('n 请输入插值区间上限 b：n'); n = input('n 输入插值节点数 n:n'); f = inline(s); x = a:(b-a)/(n-1):b; for o=1:n oo = a+(o-1)*(b-a)/(n-1); y(o) = f(oo); end end if type=2 n = input('n 输入插值节点数 n:n'); x = input('n 输入插值节点构成的向量 x_i:n'); y = input(
3、'n 输入插值节点对应的函数值构成的向量 y_i:n'); end h = zeros(n-1,1); miu = zeros(n-2,1); lambda = zeros(n-2,1); S = zeros(n-1,4); for i=1:n-1 h(i) = x(i+1)-x(i); end for j=1:n-2 miu(j) = h(j)/(h(j)+h(j+1); lambda(j) = 1-miu(j); d(j) = 6/(h(j)+h(j+1)*(y(j+2)-y(j+1)/h(j+1)-(y(j+1)-y(j)/h(j); end %边界条件的选择 bound
4、ary = input('n 选择封闭方程组的边界条件: n 第一类边界条件输入 1， n 第二类边界条件输入 2，n 第三类边界条件输入 3。n'); while boundary=1 pause(1) boundary = input('n 选择封闭方程组的边界条件: n 第一类边界条件输 入 1，n 第二类边界条件输入 2，n 第三类边界条件输入 3。n'); end %第一类边界条件 AM=D 附录 3 求三次样条插值函数的 matlab 程序 24 if boundary = 1 A = zeros(n-2,n-2); D = zeros(
5、n-2,1); M_0 = input('n 输入插值区间左端点 a 处的二阶导数值（若为自然三次 样条插值函数，输入 0） ：n'); M_n = input('n 输入插值区间右端点 b 处的二阶导数值（若为自然三次 样条插值函数，输入 0） ：n'); for k = 2:n-3 A(k,k-1:k+1) = miu(k),2,lambda(k); D(k) = d(k); end A(1,1:2) = 2,lambda(1); A(n-2,n-3:n-2) = miu(n-2),2; D(1) = d(1)-miu(1)*M_0; D(n-2) = d(
6、n-2)-lambda(n-2)*M_n; end %第二类边界条件 if boundary = 2 A = zeros(n,n); D = zeros(n,1); dydx_a = input('n 输入插值区间左端点 a 处的一阶导数值：n'); dydx_b = input('n 输入插值区间右端点 b 处的一阶导数值：n'); d_0 = 6/h(1)*(y(2)-y(1)/h(1)-dydx_a); d_n = 6/h(n-1)*(dydx_b-(y(n)-y(n-1)/h(n-1); for k = 2:n-1 A(k,k-1:k+1) = miu(
7、k-1),2,lambda(k-1); D(k) = d(k-1); end A(1,1:2) = 2,1; A(n,n-1:n) = 1,2; D(1) = d_0; D(n) = d_n; end %第三类边界条件 if boundary =3 A = zeros(n-1,n-1); D = zeros(n-1,1); miu_n = h(n-1)/(h(n-1)+h(1); lambda_n = 1-miu_n; d_n = 6/(h(n-1)+h(1)*(y(2)-y(n)/h(1)-(y(n)-y(n-1)/h(n-1); for k=2:n-2 A(k,k-1:k+1) = miu
8、(k),2,lambda(k); D(k) = d(k); end A(1,1:2) = 2,lambda(1); A(1,n-1) = miu(1); A(n-1,1) = lambda_n; A(n-1,n-2:n-1) = miu_n,2; D(1) = d(1); D(n-1) = d_n; end %用追赶法求解一、二类边界条件下的三弯矩方程组 if boundary =1 | boundary =2 N = n-(4-2*boundary); 计算方法上机报告 25 l = zeros(N,1); u = zeros(N,1); yy = zeros(N,1); M = zeros
9、(N,1); u(1) = A(1,1); yy(1) = D(1); for p=2:N l(p) = A(p,p-1)/u(p-1); u(p) = A(p,p)-l(p)*A(p-1,p); yy(p) = D(p)-l(p)*yy(p-1); end M(N) = yy(N)/u(N); q =N; while q=1 q = q-1; M(q) = (yy(q)-A(q,q+1)*M(q+1)/u(q); end if boundary =1 M = M_0;M;M_n; end end %用 LU 分解求解第三类边界条件下的三弯矩方程组 if boundary =3 l = zer
10、os(n-1,n-1); u = zeros(n-1,n-1); yy = zeros(n-1,1); u(1,:) = A(1,:); l(:,1) = A(:,1)/u(1,1); l(1,1) = 1; yy(1) = D(1); M1 = zeros(n-1,1); M = zeros(n,1); for a=2:n-2 l(a,a) = 1; lu1 = 0; for c=1:a-1 lu1 = lu1+l(a,c)*u(c,a); end u(a,a) = A(a,a)-lu1; for b=a+1:n-1 lu2 = 0; lu3 = 0; for d=1:a-1 lu2 = l
11、u2+l(a,d)*u(d,b); lu3 = lu3+l(b,d)*u(d,a); end u(a,b) = A(a,b)-lu2; l(b,a) = (A(b,a)-lu3)/u(a,a); end end lu = 0; for e=1:n-2; lu = lu+l(n-1,e)*u(e,n-1); end u(n-1,n-1) = A(n-1,n-1)-lu; l(n-1,n-1) = 1; 附录 3 求三次样条插值函数的 matlab 程序 26 for f = 2:n-1; ly = 0; for g = 1:f-1 ly = ly+l(f,g)*yy(g); end yy(f)
12、= D(f)-ly; end M1(n-1) = yy(n-1)/u(n-1,n-1); for rr=1:n-2 r = n-1-rr; uM1 = 0; for s=r+1:n-1 uM1 = uM1+u(r,s)*M1(s); end M1(r) = (yy(r)-uM1)/u(r,r); end M = M1(n-1,1);M1; end ss = 0; for t=1:n-1 S(t,1) = (M(t+1)-M(t)/(6*h(t); S(t,2) = (M(t)*x(t+1)-M(t+1)*x(t)/(2*h(t); S(t,3) = (M(t+1)*x(t)2-M(t)*x(t
13、+1)2)/(2*h(t)+(y(t+1)-y(t)/h(t)+h(t)*(M(t)-M(t+1)/6; S(t,4) = (M(t)*x(t+1)3-M(t+1)*x(t)3)/(6*h(t)+(y(t)*x(t+1)-y(t+1)*x(t)/h(t)+h(t)*(M(t+1)* x(t)-M(t)*x(t+1)/6; for x1 = x(t):(x(t+1)-x(t)/100:x(t+1) ss = ss+1; xx(ss) = x1; SS(ss) = S(t,1)*x13+S(t,2)*x12+S(t,3)*x1+S(t,4); end end plot(xx,SS,'-k&
14、#39;,'linewidth',2); hold on plot(x,y,'*k','markersize',10); hold on xlabel('x'); ylabel('S(x)'); grid; fprintf('n 所求的三次样条插值函数为：n'); for uu=1:n-1 fprintf('S(x) = %10.5f*x3+%10.5f*x2+%10.5f*x+%10.5f, %8.4f= x =%8.4fn',S(uu,1),S(uu,2),S(uu,3),S(uu,4),x(uu),x(uu+1); end
