1、 MATLAB/Simulink/C+等编程援助，请访问我的博客 编程博客：http:/ 我的QQ、Email、MSN等，点击联系方式，即可查看 联系方式：http:/ 有编程问题的朋友，请直接联系我。如果我QQ/MSN不在线， 请将问题直接发到我的邮箱 ，或者在我的博客留言，第一时间答复！ 非常荣幸能够成为大家的 QQ好友 ，博客内容不断更新中，欢迎大家持续访问！ 声明： 本资料整理于网络，仅限交流使用，切勿用做商业用途！ 浙江工业大学化材学院 jackdong A 目 录 第一章 插值方法1 1.1. Lagrange插值.1 1.2. Lagrange 插值多项式2 1.3. Newto
2、n多项式3 1.4. 切比雪夫逼近.4 1.5. 逐步插值.5 1.6. 分段三次Hermite插值6 1.7. 分段三次样条插值.7 第二章 数值积分10 2.1. 复化Simpson公式.10 2.2. 变步长梯形法.12 2.3. Romberg加速法14 2.4. 三点Gauss公式.16 第三章 常微分方程的差分解法17 3.1. 改进的Euler方法17 3.2. Heun方法18 3.3. 四次Taylor方法19 3.4. 四阶Runge-Kutta法19 3.5. Runge-Kutta-Fehlbrg法21 3.6. 二阶Adams预报校正系统.24 3.7. 改进的四阶A
3、dams预报校正系统.25 3.8. Milne-Simpson方法27 3.9. Hamming方法.29 3.10. 微分方程组四阶Runge-Kutta解法30 3.11. 线性打靶法.31 3.12. 求解三对方程组的程序.32 3.13. 有限差分法.33 附：三阶、四阶、五阶Rungkuta法.35 第四章 方程求根40 4.1. 二分法.40 4.2. 开方法.41 4.3. Newton 下山法42 4.4. 快速弦截法.43 4.5. 不动点迭代法.44 4.6. 试值法或试位法.45 4.7. Steffensen加速法46 4.8. Muller法48 浙江工业大学化材学
4、院 jackdong B 第五章 线性方程组的迭代法51 5.1. Jacobi迭代法.51 5.2. Gauss-Seidel迭代.53 5.3. 非线性Seidel迭代.56 5.4. 牛顿-拉夫森法.57 5.5. 超松弛迭代.59 5.6. 对超松弛迭代.60 第六章 线性方程组的直接法63 6.1. 追赶法.63 6.2. Cholesky方法.64 6.3. 矩阵分解方法.66 6.4. 消去法.69 第七章 数值优化71 7.1. 黄金分割法求极小值.71 7.2. 斐波那契法求极小值.73 7.3. 用2次插值求局部最小值.76 7.4. 内德-米德法求最小值.81 7.5.
5、最速下降法或梯度法.86 浙江工业大学化材学院 jackdong 1 第一章 插值方法 1.1. Lagrange插值 功能：计算多项式在x=x0处的值 - function y0,N=Lagrange_eval(X,Y,x0) % X,Y :已知的插值点坐标 % x0 :插值点 % y0 :Lagrange插值多项式在x0处的值 % N :Lagrange插值函数的权系数 m=length(X) N=zeros(m,1) y0=0; for i=1:m N(i)=1; for j=1:m if j=i N(i)=N(i)*(x0-X(j)/(X(i)-X(j); end end y0=y0+
6、Y(i)*N(i); end - 浙江工业大学化材学院 jackdong 2 1.2. Lagrange 插值多项式 功能：求基于N+1个点的拉格朗日多项式 - function C,L=lagran(X,Y) % input - X is a vector that contains a list of abscissas % - Y is a vector that contains a list of ordinates % output - C is a matrix that contains the coefficients of the lagrange interpolator
7、y polynomial - L is a matrix that contains the lagrange coefficients polynomial w=length(X); n=w-1; L=zeros(w,w); for k=1:n+1 V=1; for j=1:n+1 if k=j V=conv(V,poly(X(j)/(X(k)-X(j); end end L(k,:)=V; end C=Y*L; - 浙江工业大学化材学院 jackdong 3 1.3. Newton多项式 功能：构造和计算通过(xk,yk)=(xk,f(xk),k=0N的次数小于等于N的牛顿多项式 - fu
8、nction C,D=newpoly(X,Y) % input - X is a vector that contains a list of abscissas % - Y is a vector that contains a list of ordinates % output - C is a matrix that contains the coefficients of the newton interpolatory polynomial - D is the divided-difference table n=length(X); D=zeros(n,n); D(:,1)=Y
9、; for j=2:n for k=j:n D(k,j)=(D(k,j-1)-D(k-1,j-1)/(X(k)-X(k-j-1); end end C=D(n,n); for k=(n-1):-1:1 C=conv(c,poly(X(k); m=length(C); C(m)=C(m)+D(k.k); end - 浙江工业大学化材学院 jackdong 4 1.4. 切比雪夫逼近 功能：构造并计算-1，1上的N次切比雪夫多项式 - function C,X,Y=cheby(fun,n,a,b) % input - fun is the string function to be approxi
10、mated % -a and b are the left and right endpoints % - n is the degree of the chebyshev interpolating polynomial % output - C is the coefficient list for the polynomial - X contains the abscissas - Y contains the ordinates if nargin =2,a=1;b=1;end d=pi/(2*n+2); C=zeros(1,n+1); for k=1:n+1 X(k)=cos(2*
11、k-1)*d); end X=(b-a)*X/2+(a+b)/2; x=X; Y=eval(fun); for k=1:n+1 z=(2*k-1)*d; for j=1:n+1 C(j)=C(j)+Y(k)*cos(j-1)*z); 浙江工业大学化材学院 jackdong 5 end end C=2*C/(n+1); C(1)=C(1)/2; - 1.5. 逐步插值 功能：计算逐步插值多项式在x0处的值 - function y0=neville_eval(X,Y,x0) % X,Y :已知的插值点坐标 % x0 :插值点 % y0 :neville逐步插值多项式在x0处的值 m=length
12、(X) P=zeros(m,1) P1=zeros(m,1) P=Y; for i=1:m P1=P; k=1; for j=i+1:m k=k+1; P(j)=P(j-1)+(P1(j)-P1(j-1)*(x0-X(k-1)/(X(j)-X(k-1); 浙江工业大学化材学院 jackdong 6 end if abs(P(m)-P(m-1)=X(i) end; - 附：函数值为向量形式的复合simpson求积法 function I=simpson_n(fname,a,b,n) %调用格式：I=simpson_n(fname,a,b,n) %其中a,b为积分区间两个端点，n为子区间数目 h=
13、(b-a)/n; x=a+(0:n)*h; f=feval(fname,x); I=simpson_h(f,h) % 调用上面编译好的simpson_h函数 - 2.2. 变步长梯形法 功能：利用变步长梯形法计算函数f(x)在给定区间的积分值 - function T,n=bbct(f,a,b,eps) % f:被积函数句柄 % a,b:积分区间的两个端点 % eps:精度 % n:二分区间的次数 浙江工业大学化材学院 jackdong 13 % T:用变步长梯形法求得的积分值 h=b-a; fa=feval(f,a); fb=feval(f,b); T1=h*(fa+fb)/2; T2=T1
14、/2+h*feval(f,a+h/2)/2; n=1; %按变步长梯形法求积分值 while abs(T2-T1)=eps h=h/2; T1=T2; S=0; x=x+h/2; while xeps 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)/(4k-1); end err=abs(R(J+1,J)-R(J+1,J+1); end quad=R(J+1,J+
15、1); 浙江工业大学化材学院 jackdong 16 - 2.4. 三点Gauss公式 功能：利用三点Gauss公式计算被积函数f(x)在给定区间的积分值 - 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; = 浙江工业大学化材学院 jackdon
16、g 17 第三章 常微分方程的差分解法 3.1. 改进的Euler方法 功能：用改进Euler法求解常微分方程 - function E=MendEuler(f,a,b,n,ya) % f:微分方程右端函数句柄 % 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)=y(a); 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),
17、y1); y(i+1)=(y1+y2)/2; end E=x,y; - 浙江工业大学化材学院 jackdong 18 3.2. Heun方法 功能：通过计算yk+1=yk+(h/2)*(f(tk,yk)+f(tk+1,yk+h*f(tk,yk)求a,b上的初值问题y=f(t,y),y(a)=y0,的近似解，其中k=0,1,.M-1 - function H=henu(f,a,b,ya,M) % input - f is the function entered as a string f - a and b are the left and right end points - ya is t
18、he initial condition y(a) - M is the number of steps % output - H=T Y where T is the vector of abscissas and Y is the vector of ordinates h=(b-a)/M; T=zeros(1,M+1); Y=zeros(1,M+1); T=a:h:b; Y(1)=ya; for i=1:M k1=feval(f,T(i),Y(i); k2=feval(f,T(i+1),Y(i)+h*k1); Y(i+1)=Y(i)+(h/2)*(k1+k2); end H=T Y; -
19、 浙江工业大学化材学院 jackdong 19 3.3. 四次Taylor方法 功能：通过计算y,y,y和y，并在每一步中适用泰勒多项式，求a,b上的初值问题y=f(t,y),y(a)=y0的近似解 - function T4=taylor(df,a,b,ya,M) % input - df=y y y y entered as a string df ,where y=f(t,y) % - a and b are the left and right end points % - ya is the initial condition y(a) % - M is the number of
20、steps % output - T4=T Y where T is the vector of abscissas and Y is the vector of ordinates h=(b-a)/M; T=zeros(1,M+1); Y=zeros(1,M+1); T=a:h:b; Y(1)=ya; for i=1:M D=feval(df,T(i),Y(i); Y(i+1)=Y(i)+h*(D(1)+h*(D(2)/2+h*(D(3)/6+h*D(4)/24); end T4=T Y; - 3.4. 四阶Runge-Kutta法 浙江工业大学化材学院 jackdong 20 功能：用四阶
21、Runge-Kutta法求解常微分方程 - function R=Rungkuta4(f, a, b, n, ya) % f:微分方程右端函数句柄 % a,b:自变量取值区间的两个端点 % n:区间等分的个数 % ya:函数初值y(a) % R=x,y:自变量X和解Y所组成的矩阵 h=(b-a)/n; x=zeros(,n+1); y=zeros(1,n+1); x=a:h:b; y(1)=ya; for i=1:n k1=h*feval(f,x(i),y(i); k2=h*feval(f,x(i)+h/2,y(i)+k1/2); k3=h*feval(f,x(i)+h/2,y(i)+k2/2
22、); k4=h*feval(f,x(i)+h,y(i)+k3); y(i+1)=y(i)+(k1+2*k2+2*k3+k4)/6; end R=x,y; - 浙江工业大学化材学院 jackdong 21 3.5. Runge-Kutta-Fehlbrg法 功能：用误差控制和步长方法求解a,b上的初值问题y=f(t,y),y(a)=y0的近似解 - function R=rkf45(f,a,b,ya,M,tol) % input - f is the function entered as a string f % - a and b are the left and right end poi
23、nts % - ya is the initial condition y(a) % - M is the number of steps % -tol is the tolerance % output - R=T Y where T is the vector of abscissas and Y is the vector of ordinates a2=1/4;a3=3/8;a4=12/13;a5=1;a6=1/2; b2=1/4;b3=3/32;b4=1932/2197;b5=439/216;b6=-8/27; c3=9/32;c4=-7200/2197;c5=-8;c6=2; d4
24、=7296/2197;d5=3680/513;d6=-3544/2565; e5=-845/4104;e6=1859/4104; f6=-11/40; r1=1/360;r3=-128/4275;r4=-2197/75240;r5=1/50;r6=2/55; n1=25/216;n3=1408/2565;n4=2197/4104;n5=-1/5; big=1e15; h=(b-a)/M; hmin=h/64;hmax=64*h;max1=200; Y(1)=ya;T(1)=a; j=1; 浙江工业大学化材学院 jackdong 22 br=b-0.00001*abs(b); while (T(
25、j)br) h=b-T(j); end k1=h*feval(f,T(j),Y(j); y2=Y(j)+b2*k1; if bigbr) T(j+1)=b; else T(j+1)=T(j)+h; end j=j+1; end if (err=0) s=0; else s=0.84*(tol*h/err)0.25; end if (s2*hmin) h=h/2; end if (s1.50) else M=j; end end R=T Y; - 3.6. 二阶Adams预报校正系统 功能：用二阶Adams预报校正系统求解常微分方程 - function A=Adams(f,a,b,n,ya)
26、% f:微分方程右端函数句柄 % a,b:自变量取值区间的两个端点 % n:区间等分的个数 % ya:函数初值y(a) % A=x,y:自变量X和解Y所组成的矩阵 h=(b-a)/n; x=zeros(1,n+1); y=zeros(1,n+1); x=a:h:b; 浙江工业大学化材学院 jackdong 25 y(1)=ya; for i=1:n if i=1 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; dy1=feval(f,x(i),y(i); dy2=feval(f,x(i+1
27、),y(i+1); else y(i+1)=y(i)+h*(3*dy2-dy1)/2; P=feval(f,x(i+1),y(i+1); y(i+1)=y(i)+h*(P+dy2)/2; dy1=dy2; dy2=feval(f,x(i+1),y(i+1); end end A=x,y; - 3.7. 改进的四阶Adams预报校正系统 功能：用改进的Adams预报校正系统求解常微分方程 - function A=CAdams4PC(f,a,b,n,ya) 浙江工业大学化材学院 jackdong 26 % f:微分方程右端函数句柄 % a,b:自变量取值区间的两个端点 % n:区间等分的个数 %
28、 ya:函数初值y(a) % A=x,y:自变量X和解Y所组成的矩阵 if n4 break: end h=(b-a)/n; x=zeros(1,n+1); y=zeros(1,n+1); x=a:h:b; y(1)=ya; F=zero(1,4); for i=1:n if i4 k1=feval(f,x(i),y(i); k2=feval(f,x(i)+h/2,y(i)+(h/2)*k1); k3=feval(f,x(i)+h/2,y(i)+(h/2)*k2); k4=feval(f,x(i)+h,y(i)+h*k3); y(i+1)=y(i)+(h/6)*(k1+2*k2+2*k3+k4
29、); elseif i=4 浙江工业大学化材学院 jackdong 27 F=feval(f,x(i-3:i),y(i-3:i); py=y(i)+(h/24)*(F*-9,37,-59,55); p=feval(f,x(i+1),py); F=F(2) F(3) F(4) p; y(i+1)=y(i)+(h/24)*(F*1,5,-19,9); p=py;c=y(i+1); else F=feval(f,x(i-3:i),y(i-3:i); py=y(i)+(h/24)*(F*-9,37,-59,55); my=py-251*(p-c)/270; m=feval(f,x(i+1),my); F=F(2) F(3) F(4) m; cy=y(i)+(h/24)*(F*1,5,-19,9); y(i+1)=cy+19*(py-cy)/270; p=py;c=cy; end end A=x,y; - 3.8. Milne-Simpson方法 功能：用预报子：p(k+1)=y(k-3)+(4h/3)*(2f(k-2)-f(k-1)+2*f(k) 和校正子：y(k+1)=y(k-1)+(h/3)*(f(k-1)+4f(k)+f(k+1),求区间a,b上的初值问题
