2022年2022年计算方法上机作业求三次样条插值函数的matlab程序 .pdf

《2022年2022年计算方法上机作业求三次样条插值函数的matlab程序 .pdf》由会员分享，可在线阅读，更多相关《2022年2022年计算方法上机作业求三次样条插值函数的matlab程序 .pdf（4页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、计算方法上机报告23 附录 3 求三次样条插值函数的matlab程序clc;clear; type = input(n 选择插值节点的类型：n 插值对象为连续函数请输入1，n插值对象为列表函数请输入2。n); while type =1 & type=2 fprintf(n 输入值有误，请重新输入n); pause(1) type = input(n 选择插值节点的类型：n 插值对象为连续函数请输入1，n 插值对象为列表函数请输入2。n); end if type=1 s = input(n 请输入连续函数的表达式：nf(x) = ,s); a = input(n 请输入插值区间下限a：n);

2、 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(n 输入插值节点对应的函数值构成的向量y_i:n); end h = zeros(n-1,1); miu = zeros(n-2,1); lambda =

3、 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 %边界条件的选择boundary = input(n 选择封闭方程组的边界条件: n 第一类边界条件输入1，n 第二类边界条件输入2， n 第三类边界条件输入3。n); while boundary=1

4、& boundary=2 & boundary=3 fprintf(n 输入值有误，请重新输入n); pause(1) boundary = input(n 选择封闭方程组的边界条件: n 第一类边界条件输入 1，n 第二类边界条件输入2，n 第三类边界条件输入3。n); end %第一类边界条件AM=D 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 1 页，共 4 页 - - - - - - - - - 附录 3 求三次样条插值函数的matlab 程序24 if boundary

5、 = 1 A = zeros(n-2,n-2); D = zeros(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-

6、2) = d(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(k-1),2,lambda(k-

7、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(k),2,lambda(k);

8、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); 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 2 页，共 4

9、页 - - - - - - - - - 计算方法上机报告25 l = zeros(N,1); u = zeros(N,1); yy = zeros(N,1); M = zeros(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); e

10、nd if boundary =1 M = M_0;M;M_n; end end %用 LU 分解求解第三类边界条件下的三弯矩方程组if boundary =3 l = zeros(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

11、,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 = lu2+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; 名师资料总结 - -

12、 -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 3 页，共 4 页 - - - - - - - - - 附录 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) = 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

13、,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+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(

14、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,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名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 4 页，共 4 页 - - - - - - - - -