if x0>=X(1,j)& x0 y0=((C(4,j)*omg+C(3,j))omg+C(2,j))*omg+C(1,j); end end end 算例1：给定数据表： Xi 0.25 yi 0.5000 0.30 0.5477 0.39 0.6245 0.45 0.6708 0.53 0.7280 试求三次样条插值函数S(x)，其中：S’(0.25)=1.0000,S’(0.53)=0.6868. 解：X=[0.25,0.30,0.39,0.45,0.53];Y=[0.5000,0.5477,0.6254,0.6708,0.7280]; s0=1.0000;sN=0.6868; [y0,C]=Spline_interp(X,Y,s0,sN,x0); plot(0.25:0.01:0.30,polyval(C(:.1),0:0.01:0.05),’r-.’); hold on plot(0.30:0.01:0.39,polyval(C(:.2),0:0.01:0.09),’b’); plot(0.39:0.01:0.45,polyval(C(:.3),0:0.01:0.06),’k-*’); plot(0.45:0.01:0.53,polyval(C(:.4),0:0.01:0.08)); 第二章 数值积分 2.1 Simpson公式 利用复化Simpson公式求被积函数f(x)在给定区间上的积分值. MATLAB文件：（文件名：FSimpson.m） function S=FSimpson(f,a,b,N) %f表示被积函数句柄 %a,b表示被积区间端点 %N表示区间个数 %S是用复化Simpson公式求得的积分值 h=(b-a)/N; fa=feval(f,a); fb=feval(f,b); S=fa+fb; x=a; for i=1:N x=x+h/2; fx=feval(f,x); S=S+4*fx; x=x+h/2; fx=feval(f,x); S=S+2*fx; end S=h*S/6; 5 算例2：利用复化Simpson公式计算积分S?解：后面都要用到的f1: function f=f1(x) f=x/(4+x^2); 令f=@f1;a=0;b=1; 运行S=FSimpson(f,a,b,N) 这里的N值需要自己输入。 xdx?04?x2. 12.2 变步长梯形法 利用变步长梯形法求被积函数f(x)在给定区间上的积分值. MATLAB文件：（文件名：bbct.m） function [T,n]=bbct(f,a,b,eps) %f表示被积函数句柄 %a,b被积区间端点 %eps精度 %T是用变步长梯形法求得的积分值 %n表示二分区间的次数 h=b-a; fa=feval(f,a); fb=feval(f,b); T1=h*(fa+fb)/2; T2=T1/2+h*feval(f,a+h/2)/2; n=1; %按变步长梯形法求积分值 while abs(T2-T1)>=eps h=h/2; T1=T2; S=0; x=a+h/2; while x fx=feval(f,x); S=S+fx; x=x+h; end T2=T1/2+S*h/2; n=n+1; end T=T2; 算例3：利用变步长梯形法计算积分T?解： 6 xdx?04?x2. 1function f=f1(x) f=x/(4+x^2); 令f=@f1;a=0;b=1; 运行 [T,n]=bbct(f,a,b,eps) 这里的eps值需要自己输入。 2.3 Romberg加速算法 利用Romberg加速算法计算被积函数f(x)在给定区间上的积分值. MATLAB文件：（文件名：Romberg.m） function [quad,R]=Romberg(f,a,b,eps) %f表示被积函数句柄 %a,b被积区间端点 %eps精度 %quad是用Romberg加速算法求得的积分值 %R为Romberg表 %err表示误差的估计 h=b-a; R(1,1)=h*(feval(f,a)+feval(f,b))/2; M=1; J=0; err=1; while err>eps J=J+1; h=h/2; S=0; for p=1:M x=a+h*(2*p-1); S=S+feval(f,x); end R(J+1,1)=R(J,1)/2+h*S; M=2*M; for k=1:J R(J+1,k+1)=R(J+1,k)+(R(J+1,k)-R(J,k))/(4^k-1); end err=abs(R(J+1,J)-R(J+1,J+1)); end quad=R(J+1,J+1); 算例4：利用Romberg加速算法计算积分R?解： function f=f1(x) f=x/(4+x^2); 令f=@f1;a=0;b=1; 运行 [quad,R]=Romberg(f,a,b,eps) 这里的eps值需要自己输入。 7 xdx?04?x2. 1 2.4 三点Gauss公式 利用三点Gauss公式计算被积函数f(x)在给定区间上的积分值. MATLAB文件：（文件名：TGauss.m） function G=TGauss(f,a,b) %f表示被积函数句柄 %a,b被积区间端点 %G是用三点Gauss公式求得的积分值 x1=(a+b)/2-sqrt(3/5)*(b-a)/2; x2=(a+b)/2+sqrt(3/5)*(b-a)/2; G=(b-a)*(5*feval(f,x1)/9+8*feval(f,(a+b)/2)/9+5*feval(f,x2)/9)/2; 算例5：利用三点Gauss公式计算积分R?解： function f=f1(x) f=x/(4+x^2); 令f=@f1;a=0;b=1; 运行 G=TGauss(f,a,b) 第三章 常微分方程德差分方法 3.1 改进的Euler方法 用改进的Euler方法求解常微分方程. MATLAB文件：（文件名：MendEuler.m） function E=MendEuler(f,a,b,N,ya) %f是微分方程右端函数句柄 %a,b是自变量的取值区间[a,b]的端点 %N是区间等分的个数 %ya表初值y(a) %E=[x’,y’]是自变量X和解Y所组成的矩阵 h=(b-a)/N; y=zeros(1,N+1); x=zeros(1,N+1); y(1)=ya; x=a:h:b; for i=1:N y1=y(i)+h*feval(f,x(i),y(i)); y2=y(i)+h*feval(f,x(i+1),y1); y(i+1)=(y1+y2)/2; end E=[x’,y’]; 8 xdx?04?x2. 1
相关推荐：
loading