Matlab数值计算程序

2021-04-28 来源：爱站旅游

导读Matlab数值计算程序

浙江工业大学化材学院 jackdong

======================================================================

数值计算方法

算法设计及其 MATLAB 实现

浙江工业大学

Jackdong QQ: 7214478

Email:jackdong81413@tom.com

2006年12月17日

======================================================================

浙江工业大学化材学院 jackdong

第一章插值方法............................................................1

1.1. Lagrange插值...................................................................................................................1 1.2. Lagrange 插值多项式......................................................................................................2 1.3. Newton多项式..................................................................................................................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. 改进的四阶Adams预报校正系统.................................................................................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

浙江工业大学化材学院 jackdong

第五章线性方程组的迭代法....................................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. 最速下降法或梯度法.....................................................................................................86

浙江工业大学化材学院 jackdong

第一章插值方法

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

y0=y0+Y(i)*N(i); end

----------------------------------------------------------------------------------------------------------------------

浙江工业大学化材学院 jackdong

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 interpolatory 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

L(k,:)=V; end

C=Y*L;

----------------------------------------------------------------------------------------------------------------------

浙江工业大学化材学院 jackdong

1.3. Newton多项式

功能：构造和计算通过(xk,yk)=(xk,f(xk)),k=0....N的次数小于等于N的牛顿多项式 ----------------------------------------- function [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';

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

1.4. 切比雪夫逼近

功能：构造并计算[-1，1]上的N次切比雪夫多项式 ----------------------------------------- function [C,X,Y]=cheby(fun,n,a,b)

% input - fun is the string function to be approximated

% -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*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

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(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

end

if abs(P(m)-P(m-1))<10^-6

y0=P(m);

return;

end end

y0=P(m);

----------------------------------------------------------------------------------------------------------------------

1.6. 分段三次Hermite插值

功能：利于分段三次Hermite插值计算插值点处的函数近似值 -----------------------------------------

function y0=Hermite_inter(X,Y,DY,x0)

% X,Y :已知的插值点坐标向量

% x0 :插值点

% y0 :Hermite插值多项式在x0处的值

% DY :插值点处的导数值

N＝length(X)

for i=1:N

if x0 >=X(i) & x0 <=X(i+1)

k=i;break;

end end

a1=x0-X(k+1);

浙江工业大学化材学院 jackdong

a2=x0-X(k);

a3=X(k)-X(k+1);

y0=(a1/a3)^2*(1-2*a2/a3)*Y(k)+(-a2/a3)^2*(1+2*a1/a3)*Y(k+1)+(a1/a3)^2*a2*DY(k)+(-a2/a3)^2*a1*DY(k+1);

----------------------------------------------------------------------------------------------------------------------

1.7. 分段三次样条插值

功能：计算在插值点处的函数值，并拟合曲线 -----------------------------------------

function [y0,C]=Spline_interp(X,Y,s0,sN,x0)

% X,Y :已知的插值点坐标向量

% s0,sN:两个端点的一次导数值

% x0 :插值点

% y0 :三次样条函数在x0处的值 % C :分段三次样条函数的系数

N=length(X)

C=zeros(4,N-1); h=zeros(1,N-1); mu=zeros(1,N-1); lmt=zeros(1,N-1);d=zeros(1,N);

h=X(1,2:N)-X(1,1:N-1);

mu(1,N-1)=1;

lmt(1,1)=1;

mu(1,2:N-2)=h(1,1:N-2)/(h(1,1:N-2)+h(1,1:N-1));

lmt(1,2:N-1)=h(1,2:N-1)/(h(1,1:N-2)+h(1,2:N-1));

d(1,1)=6*(Y(1,2)-Y(1,1))/h(1,1)-s0/h(1,1);

d(1,N)=6*(sN-(Y(1,N)-Y(1,N-1))/h(1,N-1)/h(1,N-1);

浙江工业大学化材学院 jackdong

d(1,2:N-1)=6*((Y(1,3:N)-Y(1,2:N-1))/h(1,2:N-1)-(Y(1,2:N-1)-Y(1,1:N-2))/h(1,1:N-2)/h(1,1:N-2)+h(1,2:N-1);

%追赶法解三对角方程组

bit=zeros(1,N-1);

bit(1,1)=lmt(1,1)/2;

for i=2:N-1

bit(1,i)=lmt(1,i)/(2-mu(1,i-1)*bit(1,i-1)); end

y=zeros(1,N);

y(1,1)=d(1,1)/2;

for i=2:N

y(1,i)=(d(1,i)-mu(1,i-1)*y(1,i-1))/(2-mu(1,i-1)*bit(1,i-1)); end

x=zeros(1,N);

x(1,N)=y(1,N);

for i=N-1:-1:1

x(1,i)=y(1,i)-bit(1,i)*x(1,i+1); end

v=zeros(1,N-1);

v(1,1:N-1)=(Y(1,2:N)-Y(1,1:N-1))/h(1,1:N-1);

C(4,:)=Y(1,1:N-1);

C(3,:)=v-h*(2*x(1,N-1)+x(1,2:N))/6;

浙江工业大学化材学院 jackdong

C(2,:)=x(1,1:N-1)/2;

C(1,:)=(x(1,2:N)-x(1,1:N-1))/(6*h);

if nargin<5

y0=0; else

for j=1:N-1

if x0>=X(1,j) & x0omg=x0-X(1,j);

y0=(C(4,:)*omg+C(3,:)*omg+C(2,:)*omg)+C(1,j);

end end end

======================================================================

浙江工业大学化材学院 jackdong

第二章数值积分

2.1. 复化Simpson公式

功能：利用复化Simpson公式计算被积函数f(x)在给定区间上的积分值 ----------------------------------------- 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

浙江工业大学化材学院 jackdong

S=h*S/6;

----------------------------------------------------------------------------------------------------------------------

附：函数值为向量形式的simpson求积法

function I=simpson_h(f,h)

%调用格式I＝simpson(f,h)

%f为一向量，指定已知节点处的函数值

%h为步长

n=length(f)-1;

if n==1

fprintf('Data has only one interval'),return; end;

if n==2

I=(h/3)*(f(1)+4*f(2)+f(3));return; end;

if n==3

I=(3*h/8)*(f(1)+3*f(2)+3*f(3)+f(4));return; end; I=0; if 2*floor(n)~=n % floor is a function round towards -inf

I=3*(h/8)*(f(n-2)+3*f(n-1)+3*f(n)+f(n+1));

m=n-3;

else

m=n;

浙江工业大学化材学院 jackdong

end;

I=I+(h/3)*(f(1)+4*sum(f(2:2:m))+f(m+1));

if m>2

I=I+(h/3)*2*sum(f(3:2:m)); end;

---------------------------------------------------------------------------------------------------------------------- 附：函数值为向量形式的复合simpson求积法

function I=simpson_n(fname,a,b,n)

%调用格式：I=simpson_n('fname',a,b,n)

%其中a,b为积分区间两个端点，n为子区间数目

h=(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

% T:用变步长梯形法求得的积分值

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=x+h/2;

while xfx=feval(f,x);

S=S+fx;

x=x+h;

end

T2=T1/2+S*h/2;

n=n+1; end

T=T2;

浙江工业大学化材学院 jackdong

---------------------------------------------------------------------------------------------------------------------- 附：函数值为向量形式的梯形求积法

function I=trapez_h(f,h)

%调用格式I=trapez_h(f,h)

%其中f为一向量，即给定的一系列纵坐标值

%h为步长

I=h*(sum(f)-(f(1)+f(length(f)))/2);

---------------------------------------------------------------------------------------------------------------------- 附：函数值为向量形式的复合梯形求积法

function I=trapez_n(fname,a,b,n)

%调用格式 trapez_n('fname',a,b,n)

%a，b为区间的两个端点，n为子区间个数

%其中需要调用另外一个函数：trapez_h,是梯形求积法的函数，如上。

h=(b-a)/n;

x=a+(0:n)*h;

f=feval(fname,x);

I=trapez_h(f,h)

----------------------------------------------------------------------------------------------------------------------

2.3. Romberg加速法

功能：利用Romberg加速算法计算被积函数f(x)在给定区间的积分值 ----------------------------------------- funvtion [quad,R]=Romberg(f,a,b,eps)

% f:被积函数句柄

浙江工业大学化材学院 jackdong

% a,b:积分区间的两个端点

% eps:精度

% quad:用Romberg加速算法求得的积分值

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);

浙江工业大学化材学院 jackdong

----------------------------------------------------------------------------------------------------------------------

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;

======================================================================

浙江工业大学化材学院 jackdong

第三章常微分方程的差分解法

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),y1);

y(i+1)=(y1+y2)/2; end

E=[x',y'];

----------------------------------------------------------------------------------------------------------------------

浙江工业大学化材学院 jackdong

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 the 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'];

----------------------------------------------------------------------------------------------------------------------

浙江工业大学化材学院 jackdong

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 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

功能：用四阶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);

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

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 points

% - 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=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

br=b-0.00001*abs(b);

while (T(j)if ((T(j)+h)>br)

h=b-T(j);

end

k1=h*feval(f,T(j),Y(j));

y2=Y(j)+b2*k1;

if bigk2=h*feval(f,T(j)+a2*h,y2);

y3=Y(j)+b3*k1+c3*k2;

if bigk3=h*feval(f,T(j)+a3*h,y3);

y4=Y(j)+b4*k1+c4*k2+d4*k3;

if bigk4=h*feval(f,T(j)+a4*h,y4);

y5=Y(j)+b5*k1+c5*k2+d5*k3+e5*k4;

if bigk5=h*feval(f,T(j)+a5*h,y5);

y6=Y(j)+b6*k1+c6*k2+d6*k3+e6*k4+f6*k5;

if bigk6=h*feval(f,T(j)+a6*h,y6);

err=abs(r1*k1+r3*k3+r4*k4+r5*k5+r6*k6);

浙江工业大学化材学院 jackdong

ynew=Y(j)+n1*k1+n3*k3+n4*k4+n5*k5;

if ((errY(j+1)=ynew;

if ((T(j)+h)>br)

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 ((s<0.75) & (h>2*hmin))

h=h/2;

end

if ((s>1.50) & (2*hh=2*h;

end

if ((big23

浙江工业大学化材学院 jackdong

M=j;

if (b>T(j))

M=j+1;

else

M=j;

end end

R=[T' Y'];

----------------------------------------------------------------------------------------------------------------------

3.6. 二阶Adams预报校正系统

功能：用二阶Adams预报校正系统求解常微分方程 ---------------------------------------------------- function A=Adams(f,a,b,n,ya)

% 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

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),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

% f:微分方程右端函数句柄

% a,b:自变量取值区间的两个端点

% n:区间等分的个数

% ya:函数初值y(a)

% A=[x',y']:自变量X和解Y所组成的矩阵

if n<4

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 i<4

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);

elseif i==4

浙江工业大学化材学院 jackdong

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]上的初值问题

浙江工业大学化材学院 jackdong

y'=f(t,y),y(a)=y0的近似解

-------------------------------------------------- function M=miline(f,T,Y)

% input - f is the function entered as a string 'f'

% - T is the vector of abscissas

% - Y is the vector of ordinates

% Remark: The first four coordinates of T and Y must have starting values obtained with RK4

% output - M=[T' Y'] where T is the vector of abscissas and Y is the vector of ordinates

n=length(T);

if n<5,break;end

F=zeros(1,4);

F=feval(f,T(1:4),Y(1:4));

h=T(2)-T(1);

pold=0;

yold=0;

for k=4:n-1

pnew=Y(k-3)+(4*h/3)*(F(2:4)*[2 -1 2]');

pmod=pnew+28*(yold-pold)/29;

T(k+1)=T(1)+h*k;

F=[F(2) F(3) F(4) feval(f,T(k+1),pmod)];

Y(k+1)=Y(k-1)+(h/3)*(F(2:4)*[1 4 1]');

pold=pnew;

yold=Y(k+1);

浙江工业大学化材学院 jackdong

F(4)=feval(f,T(k+1),Y(k+1)); end

M=[T' Y'];

----------------------------------------------------------------------------------------------------------------------

3.9. Hamming方法

功能：用预报子：p(k+1)=y(k-3)+(4h/3)*(2f(k-2)-f(k-1)+2*f(k))和校正子：y(k+1)=(y(k)-y(k-2))/8+(3h/8)*(-f(k-1)+2f(k)+f(k+1)),求区间[a,b]上的初值问题y'=f(t,y),y(a)=y0的近似解

---------------------------------------------------------------- function H=hamming(f,T,Y)

% input - f is the function entered as a string 'f'

% - T is the vector of abscissas

% - Y is the vector of ordinates

% Remark: The first four coordinates of T and Y must have starting values obtained with RK4

% output - H=[T' Y'] where T is the vector of abscissas and Y is the vector of ordinates

n=length(T);

if n<5,break;end

F=zeros(1,4);

F=feval(f,T(1:4),Y(1:4));

h=T(2)-T(1);

pold=0;

cold=0;

for k=4:n-1

浙江工业大学化材学院 jackdong

pnew=Y(k-3)+(4*h/3)*(F(2:4)*[2 -1 2]');

pmod=pnew+112*(cold-pold)/121;

T(k+1)=T(1)+h*k;

F=[F(2) F(3) F(4) feval(f,T(k+1),pmod)];

cnew=(9*Y(k)-Y(k-2)+3*h*(F(2:4)*[-1 2 1]'))/8;

Y(k+1)=cnew+9*(pnew-cnew)/121;

pold=pnew;

cold=cnew;

F(4)=feval(f,T(k+1),Y(k+1)); end

H=[T' Y'];

----------------------------------------------------------------------------------------------------------------------

3.10. 微分方程组四阶Runge-Kutta解法

功能：求区间[a,b]上的微分方程组

x1'(t)=f1(t,x1(t),......xn(t))

. . . . . . . . . xn'(t)=fn(t,x1(t),......xn(t))

--------------------------------------------- function [T,Z]=rks4(F,a,b,Za,M)

% input - F is the function entered as a string 'f'

% - a and b are the left and right end points

% - Za=[x(a) y(a)] are the initial conditions

浙江工业大学化材学院 jackdong

% - M is the number of steps

% output - T is the vector of steps

- Z=[x1(t) ......xn(t)] where xk(t) is the approximation to the kth dependent variable

h=(b-a)/M;

T=zeros(1,M+1);

Z=zeros(M+1,length(Za));

T=a:h:b;

Z(1,:)=Za;

for i=1:M

k1=h*feval(F,T(i),Z(i));

k2=h*feval(F,T(i)+h/2,Z(i,:)+k1/2);

k3=h*feval(F,T(i)+h/2,Z(i,:)+k2/2);

k4=h*feval(F,T(i)+h,Z(i,:)+k3);

Z(i+1,:)=Z(i,:)+(k1+2*k2+2*k3+k4)/6; end

----------------------------------------------------------------------------------------------------------------------

3.11. 线性打靶法

功能: 求区间[a,b]上的边值问题x''=p(t)*x'(t)+q(t)*x(t)+r(t),x(a)=alpha,x(b)=beta的近似解，调用四阶RK方法

-------------------------------------------------- function L=linsht(F1,F2,a,b,alpha,beta,M)

% input - F1 and F2 are the system of first-order equations representing the I.V.P,'s (9)and (10),respectively,input as

浙江工业大学化材学院 jackdong

% strings 'F1' 'F2'

% - a and b are the left and right end points

% - alpha=x(a) and beta=x(b) :boundary conditions

% - M is the number of steps

% output - L=[T' x] where T' is the (M+1)*1 vector of ordinates

%solve the system F1

Za=[alpha,0];

[T,Z]=rks4(F1,a,b,Za,M);

U=Z(:,1);

%solve the system F2

Za=[0,1];

[T,Z]=rks4(F2,a,b,Za,M);

V=Z(:,1);

%calculate the solution to the boundary value problem

X=U+(beta-U(M+1)*V)/V(M+1);

L=[T' X];

----------------------------------------------------------------------------------------------------------------------

3.12. 求解三对方程组的程序

功能：求解三对角方程组CX=B,其中C为三对角矩阵

--------------------------------------------------- function X=trisys(A,D,C,B)

%Input - A is the subdiagnonal of the coefficient matrix

浙江工业大学化材学院 jackdong

% - D is the main diagnonal of the coefficient matrix

% - C is the superdiagnonal of the coefficient matrix

% - B is the constant vector of the linear system

%output - X is the solution vector

N=length(B);

for k=2:N

mult=A(k-1)/D(k-1);

D(k)=D(k)-mult*C(k-1);

B(k)=B(k)-mult*B(k-1); end

X(N)=B(N)/D(N);

for k=N-1:-1:1

X(k)=(B(k)-C(k)*X(k+1))/D(k); end

----------------------------------------------------------------------------------------------------------------------

3.13. 有限差分法

功能：用O(h^2)阶的有限差分方法求区间[a，b]上的边值问题x''=p(t)*x'(t)+q(t)*x(t)+r(t),x(a)=alpha,x(b)=beta的逼近。 -------------------------------------------------------- function F=linsht(p,q,r,a,b,alpha,beta,N)

% input - p,q and r are the coefficient function,input as strings 'p' ,'q' , 'r'

% - a and b are the left and right end points

% - alpha=x(a) and beta=x(b) :boundary conditions

% - N is the number of steps

浙江工业大学化材学院 jackdong

% output - F=[T' X'] where T' is the 1*N vector of abscissas and X'is the 1*N vector of ordinates

T=zeros(1,N+1);

X=zeros(1,N+1);

Va=zeros(1,N-2);

Vb=zeros(1,N-1);

Vc=zeros(1,N-2);

Vd=zeros(1,N-1);

h=(b-a)/N;

Vt=a+h:a+h(N-1);

Vb=-h^2*feval(r,Vt);

Vb(1)=Vb(1)+(1+h/2*feval(p,Vt(1)))*alpha;

Vb(N-1)=Vb(N-1)+(1-h/2*feval(p,Vt(N-1)))*beta;

Vd=2+h^2*feval(q,Vt);

Vta=Vt(1,2:N-1);

Va=-1-h/2*feval(p,Vta);

Vtc=Vt(1,1:N-2);

Vc=-1+h/2*feval(p,Vtc);

X=trisys(Va,Vd,Vc,Vb);

T=[a,Vt,b];

X=[alpha,X,beta];

F=[T' X];

浙江工业大学化材学院 jackdong

附：三阶、四阶、五阶Rungkuta法

三阶Rungkuta法

function Y = Rungkuta3(odefun,tspan,y0,varargin)

% Rungkuta3 Solve differential equations with a non-adaptive method of order 3.

% Y = Rungkuta3 (ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN] integrates % the system of differential equations y' = f(t,y) by stepping from T0 to % T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector. % The vector Y0 is the initial conditions at T0. Each row in the solution % array Y corresponds to a time specified in TSPAN.

% Y = Rungkuta3 (ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters % P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...).

% This is a non-adaptive solver. The step sequence is determined by TSPAN % but the derivative function ODEFUN is evaluated multiple times per step.

% The solver implements the Bogacki-Shampine Runge-Kutta method of order 3. % Example

% tspan = 0:0.1:20;

% y = Rungkuta3 (@vdp1,tspan,[2 0]); % plot (tspan,y(:,1));

% solves the system y' = vdp1 (t,y) with a constant step size of 0.1, % and plots the first component of the solution.

if ~isnumeric(tspan)

error('TSPAN should be a vector of integration steps.'); end

if ~isnumeric(y0)

error ('Y0 should be a vector of initial conditions.'); end

h = diff (tspan);

if any(sign(h(1))*h <= 0)

error('Entries of TSPAN are not in order.') end try

f0 = feval(odefun,tspan(1),y0,varargin{:}); catch

msg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr]; error(msg); end

浙江工业大学化材学院 jackdong

y0 = y0(:); % Make a column vector. if ~isequal(size(y0),size(f0))

error('Inconsistent sizes of Y0 and f(t0,y0).'); end

neq = length(y0); N = length(tspan); Y = zeros(neq,N); F = zeros(neq,3);

Y(:,1) = y0; for i = 2:N ti = tspan(i-1); hi = h(i-1); yi = Y(:,i-1);

F(:,1) = feval(odefun,ti,yi,varargin{:});

F(:,2) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,1),varargin{:}); F(:,3) = feval(odefun,ti+0.75*hi,yi+0.75*hi*F(:,2),varargin{:}); Y(:,i) = yi + (hi/9)*(2*F(:,1) + 3*F(:,2) + 4*F(:,3)); end Y = Y.';

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 四阶龙格－库塔法

function Y =Rungkuta4(odefun,tspan,y0,varargin)

% Rungkuta4 Solve differential equations with a non-adaptive method of order 4.

% Y = Rungkuta4 (ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN] integrates % the system of differential equations y' = f(t,y) by stepping from T0 to % T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector. % The vector Y0 is the initial conditions at T0. Each row in the solution % array Y corresponds to a time specified in TSPAN.

% Y = Rungkuta4 (ODEFUN, TSPAN, Y0, P1, P2...) passes the additional parameters % P1, P2... to the derivative function as ODEFUN (T, Y, P1, P2...).

% This is a non-adaptive solver. The step sequence is determined by TSPAN % but the derivative function ODEFUN is evaluated multiple times per step. % The solver implements the classical Runge-Kutta method of order 4. % Example

% tspan = 0:0.1:20;

% y = Rungkuta4 (@vdp1,tspan,[2 0]); % plot(tspan,y(:,1));

% solves the system y' = vdp1(t,y) with a constant step size of 0.1, % and plots the first component of the solution.

浙江工业大学化材学院 jackdong

if ~isnumeric(tspan)

error('TSPAN should be a vector of integration steps.'); end

if ~isnumeric(y0)

error('Y0 should be a vector of initial conditions.'); end

h = diff(tspan);

if any(sign(h(1))*h <= 0)

error('Entries of TSPAN are not in order.') end try

f0 = feval(odefun,tspan(1),y0,varargin{:}); catch

msg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr]; error(msg); end

y0 = y0(:); % Make a column vector. if ~isequal(size(y0),size(f0))

error('Inconsistent sizes of Y0 and f(t0,y0).'); end

neq = length(y0); N = length(tspan); Y = zeros(neq,N); F = zeros(neq,4);

Y(:,1) = y0; for i = 2:N ti = tspan(i-1); hi = h(i-1); yi = Y(:,i-1);

F(:,1) = feval(odefun,ti,yi,varargin{:});

F(:,2) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,1),varargin{:}); F(:,3) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,2),varargin{:}); F(:,4) = feval(odefun,tspan(i),yi+hi*F(:,3),varargin{:}); Y(:,i) = yi + (hi/6)*(F(:,1) + 2*F(:,2) + 2*F(:,3) + F(:,4)); end Y = Y.';

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

浙江工业大学化材学院 jackdong

五阶龙格－库塔法

function Y = Rungkuta5(odefun,tspan,y0,varargin)

% Rungkuta5Solve differential equations with a non-adaptive method of order 5.

% Y = Rungkuta5 (ODEFUN, TSPAN, Y0) with TSPAN = [T1, T2, T3 ... TN] integrates % the system of differential equations y' = f(t,y) by stepping from T0 to % T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector. % The vector Y0 is the initial conditions at T0. Each row in the solution % array Y corresponds to a time specified in TSPAN.

% Y = Rungkuta5 (ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters % P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...).

% This is a non-adaptive solver. The step sequence is determined by TSPAN % but the derivative function ODEFUN is evaluated multiple times per step. % The solver implements the Dormand-Prince method of order 5 in a general % framework of explicit Runge-Kutta methods. % Example

% tspan = 0:0.1:20;

% y = Rungkuta5 (@vdp1,tspan,[2 0]); % plot(tspan,y(:,1));

% solves the system y' = vdp1(t,y) with a constant step size of 0.1, % and plots the first component of the solution.

if ~isnumeric(tspan)

error('TSPAN should be a vector of integration steps.'); end

if ~isnumeric(y0)

error('Y0 should be a vector of initial conditions.'); end

h = diff(tspan);

if any(sign(h(1))*h <= 0)

error('Entries of TSPAN are not in order.') end try

f0 = feval(odefun,tspan(1),y0,varargin{:}); catch

msg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr]; error(msg); end

y0 = y0 (:); % Make a column vector.

浙江工业大学化材学院 jackdong

if ~isequal(size(y0),size(f0))

error('Inconsistent sizes of Y0 and f(t0,y0).'); end

neq = length(y0); N = length (tspan); Y = zeros (neq,N);

% Method coefficients -- Butcher's tableau % C | A % --+--- % | B

C = [1/5; 3/10; 4/5; 8/9; 1];

A = [ 1/5, 0, 0, 0, 0;3/40, 9/40, 0, 0, 0;44/45 -56/15, 32/9, 0, 0;... 19372/6561, -25360/2187, 64448/6561, -212/729, 0; 9017/3168, -355/33, 46732/5247, 49/176, -5103/18656];

B = [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84];

% More convenient storage A = A.'; B = B(:);

nstages = length(B); F = zeros (neq,nstages);

Y(:,1) = y0; for i = 2:N ti = tspan(i-1); hi = h(i-1); yi = Y(:,i-1);

% General explicit Runge-Kutta framework F(:,1) = feval(odefun,ti,yi,varargin{:}); for stage = 2:nstages

tstage = ti + C(stage-1)*hi;

ystage = yi + F(:,1:stage-1)*(hi*A(1:stage-1,stage-1)); F(:,stage) = feval(odefun,tstage,ystage,varargin{:}); end

Y(:,i) = yi + F*(hi*B); end Y = Y.';

======================================================================

浙江工业大学化材学院 jackdong

第四章方程求根

4.1. 二分法

功能：用二分法求解非线性方程f(x)=0在区间[a,b]内的根

------------------------------------------------------- function [x,k]=demimethod(a,b,f,emg)

% a,b :求解区间的两个端点

% f :所求方程的函数名

% emg ：精度指标

% x：所求近似解

% k: 循环次数

fa=feval(f,a);

fab=feval(f,(a+b)/2); k=0;

while abs(b-a)>emg

if fab==0

x=(a+b)/2;

return;

elsif fa*fab<0

b=(a+b)/2;

else

a=(a+b)/2;

浙江工业大学化材学院 jackdong

end

fa=feval(f,a);

fab=feval(f,(a+b)/2);

k=k+1; end

x=(a+b)/2;

----------------------------------------------------------------------------------------------------------------------

4.2. 开方法

功能：求实数的开方运算 ----------------------------

function y=Kaifang(a,eps,x0)

% a：被开放数

% eps :精度指标

% x0 :初值

% y:a的开放

x(1)=x0;

x(2)=(x(1)+a/x(1))/2; k=2;

while abs(x(k)-x(k-1))>eps

x(k+1)=(x(k)+a/x(k))/2;

k=k+1; end

y=x';

浙江工业大学化材学院 jackdong

----------------------------------------------------------------------------------------------------------------------

4.3. Newton 下山法

功能：用牛顿下山法求解非线性方程f(x)=0的根 ----------------------------------------------- function [x,k]=Mendnewton(f,x0,emg)

% f ：非线性方程

% x0 :迭代初值

% emg :精度指标

% k,u : 迭代次数和下山因子

[f1,d1]=feval(f,x0); k=1;

x(1)=x0;

x(2)=x(1)-f1/d1;

while abs(f1)>emg

u=1;

k=k+1;

[f1,d1]=feval(f,x(k));

x(k+1)=x(k)-u*f1/d1;

[f2,d2]=feval(f,x(k+1));

while abs(f2)>abs(f1)

u=u/2;

x(k+1)=x(k)-u*f1/d1;

浙江工业大学化材学院 jackdong

[f2,d2]=feval(f,x(k+1));

end end

------------------------------------------------------------------------------------------------------------------ 4.4. 快速弦截法

功能：有快速弦截法求解非线性方程f(x)=0的根 ---------------------------------------------------- function [x,k]=Fast_chord(f,x1,x2,emg)

% f:非线性方程函数

% x1,x2 : 迭代初值

% emg :精度指标

% k :循环次数 k=1;

y1=feval(f,x1);

y2=feval(f,x2);

x(k)=x2-(x2-x1)*y2/(y2-y1);

y(k)=feval(f,x(k));

k=k+1;

x(k)=x(k-1)-(x(k-1)-x2)*y(k-1)/(y(k-1)-y2);

while abs(x(k)-x(k-1))>emg

y(k)=feval(f,x(k));

x(k+1)=x(k)-(x(k)-x(k-1))*y(k)/(y(k)-y(k-1));

k=k+1; end

----------------------------------------------------------------------------------------------------------------------

浙江工业大学化材学院 jackdong

4.5. 不动点迭代法

功能：求方程x=g(x)的近似值，起始值为p0，迭代式：p(n+1)=g(p(n))

---------------------------------------------------------------------------------------------------------------------- function [k,p,err,P]=fixpt(g,p0,tol,max1)

% input - g is the iteration function input as a string 'g'

% - p0 is the initial guess for the fixed point

% - tol is the tolerance

% - max1 is the maximum number of interations

%outpt - k is the number of iterations that were carried out

% - p is the approximation to the fixed point

% - err is the error in the approximation

% - P contains the sequence {pn}

P(1)=p0;

for k=2:max1

P(k)=feval(g,P(k-1));

err=abs(P(k)-P(k-1));

relerr=err/(abs(P(k))+eps);

p=P(k);

if (errif k==max1

disp('maximum number of iterations exceeded') end P=P';

----------------------------------------------------------------------------------------------------------------------

浙江工业大学化材学院 jackdong

4.6. 试值法或试位法

功能：求解方程f(x)=0在[a,b]内的根，前提是f(x)连续，且f(a)*f(b)<0

-------------------------------------------------------

function [c,err,yc]=regula(f,a,b,delta,epsilon,max1)

% input - f is the iteration function input as a string 'f'

% - a and b are the left and right points

% - delta is the tolerance for the zero

% - epsilon is the tolerance for the value of f at the zero

% - max1 is the maximum number of interations

%outpt - c is the zero

% - err is the error estimate for c

% - yc=f(c)

ya=feval(f,a);

yb=feval(f,b);

if ya*yb>0

disp('Note f(a)*f(b)>0'),

break; end

for k=1:max1

dx=yb*(b-a)/(yb-a);

c=b-dx;

ac=c-a;

yc=feval(f,c);

浙江工业大学化材学院 jackdong

if yc==0 ,break;

elseif yb*yc<0

b=c; else

a=c; ya=yc; end

dx=min(abs(dx),ac);

if abs(dx)if abs(yc)err=abs(b-a)/2; yc=feval(f,c);

---------------------------------------------------------------------------------------------------------------------

4.7. Steffensen加速法

功能：给定初始近似值p0，快速寻找不动点方程x=g(x)的解，假设g(x)和g'(x)是连续的，且|'g(x)|<1,通常的不动点迭代缓慢线性收敛到p ------------------------------------------

function [P,Q]=steff(f,df,p0,delta,epsilon,max1)

% input - f is the object function input as a string 'f'

% - df is the derivative of f input as a string 'f'

% - delta is the tolerance for p0

% - epsilon is the tolerance for the function value of y

% - max1 is the maximum number of interations

%outpt - P is the Steffensen approximation to the zero

% - Q is the matrix containing the Steffensen sequence

浙江工业大学化材学院 jackdong

R=zeros(max1,3);

R(1,1)=p0;

for k=1:max1

for j=2:3

nrdenom=feval(df,R(k,j-1));

if nrdenom==0

disp('division by zero in Newton method');

break;

else

R(k,j)=R(k,j-1)-feval(f,R(k,j-1))/nrdenom;

end

aadenom=R(k,3)-2*R(k,2)+R(k,1);

if aadenom==0

disp('division by zero in Aitken's Acceleration);

break;

else

R(k+1,)=R(k,1)-(R(k,2)-R(k,1))^2/aadenom;

end

if (nrdenom==0) | (aadenom==0)

break;

end

浙江工业大学化材学院 jackdong

err=abs(R(k,1)-R(k+1,1));

relerr=err/(abs(R(k+1,1))+delta);

y=feval(f,R(k+1,1));

if (errP=R(k+1,1);

Q=R(1:k+1,:);

break;

end end

----------------------------------------------------------------------------------------------------------------------

4.8. Muller法

功能：给定3个初始近似值p0，p1,p2，求方程f(x)=0的根 ----------------------------------------------------------------- function [p,y,err]=muller(f,p0,p1,p2,delta,epsilon,max1)

% input - f is the object function input as a string 'f'

% - p0,p1,p2 are the initial approximations

% - delta is the tolerance for p0,p1,p2

% - epsilon is the tolerance for the function values y

% - max1 is the maximum number of interations

%outpt - p is the muller approximation to the zero of f

% - y is the function value y=f(p)

% - err is the error in the approximation of p

P=[p0 p1 p2];

浙江工业大学化材学院 jackdong

Y=feval(f,P);

for k=1:max1

h0=P(1)-P(3);

h1=P(2)-P(3);

e0=Y(1)-Y(3);

e1=Y(2)-Y(3);

c=Y(3);

denom=h1*h0^2-h0*h1^2;

a=(e0*h1-e1*h0)/denom;

b=(e1*h0^2-e0*h1^2)/denom;

if b^2-4*a*c>0

disc=sqrt(b^2-4*a*c);

else

disc=0;

end

if b<0

disc=-disc;

end

z=-2*c/(b+disc);

p=P(3)+z;

if abs(p-P(2))Q=[P(2) P(1) P(3)];

浙江工业大学化材学院 jackdong

P=Q;

Y=feval(f,P);

end

if abs(p-P(3))R=[P(1) P(3) P(1)];

P=R;

Y=feval(f,P);

end

P(3)=P;

Y(3)=feval(f,P(3));

y=Y(3);

err=abs(z);

relerr=err/(abs(p)+dela);

if (errbreak;

end end

浙江工业大学化材学院 jackdong

第五章线性方程组的迭代法

5.1. Jacobi迭代法

功能：用Jacobi迭代法求解线性方程组A*x=b

-----------------------------------------------

function [x,k]=Jacobimethod(A,b,x0,N,emg)

% A:线性方程组左端矩阵

% b:线性方程组右端向量

% x0:迭代初值

% N:迭代次数上界，若迭代次数大于n，则迭代失败

% emg:精度指标

% k:迭代次数

% x:用迭代法求得的线性方程组的近似解

n=length(A);

x1=zeros(n,1);x2=zeros(n,1);

x1=x0; k=0;

r=max(abs(b-A*x1));

while r>emg

for i=1:n

sum=0;

for j=1:n

if i~=j

浙江工业大学化材学院 jackdong

sum=sum+A(i,j)*x1(j);

end

x2(i)=(b(i)-sum)/A(i,i);

end

r=max(abs(x2-x1));

x1=x2;

k=k+1;

if k>N

disp('迭代失败，返回')；

return;

end end

----------------------------------------------------

附：Jacobi迭代程序的另一种写法

function X=jacobi(A,B,P,delta,max1)

% input - A is an N*N nonsingular matrix

- B is an N*1 matrix

- P is an N*1 matrix :the initial guess

- delta is the tolerance for P

- max1 is the maximum number of iterations

% output - X is an N*1 matrix : the jacobi approximation to the solution of AX=B

浙江工业大学化材学院 jackdong

N=length(B);

for k=1:max1

for j=1:N

X(j)=(B(j)-A(j,[1;j-1,j+1:N])*P([1:j+1,j+1:N]))/A(j,j);

end

err=abs(norm(X'-P));

relerr=err/(norm(X)+eps);

P=X';

if (errbreak;

end end

X=X';

----------------------------------------------------------------------------------------------------------------------

5.2. Gauss-Seidel迭代

功能：用Gauss-Seidel迭代求解线性方程组A*x=b

---------------------------------------------------------------------------------------------------------------------- function [x,k]=Gaussmethod(A,b,x0,N,emg)

% A:线性方程组左端矩阵

% b:线性方程组右端向量

% x0:迭代初值

% N:迭代次数上界，若迭代次数大于n，则迭代失败

% emg:精度指标

浙江工业大学化材学院 jackdong

% k:迭代次数

% x:用迭代法求得的线性方程组的近似解

n=length(A);

x1=zeros(n,1);x2=zeros(n,1);

x1=x0;

r=max(abs(b-A*x1)); k=0;

while r>emg

for i=1:n

sum=0;

for j=1:n

if j>i

sum=sum+A(i,j)*x1(j);

elseif jsum=sum+A(i,j)*x2(j);

end

x2(i)=(b(i)-sum)/A(i,i);

end

r=max(abs(x2-x1));

x1=x2;

k=k+1;

浙江工业大学化材学院 jackdong

if k>N

disp('迭代失败，返回')；

return;

end end

x=x1;

---------------------------------------------------------------- 附：Gauss-Seidel迭代程序的另一种写法

function X=gseid(A,B,P,delta,max1)

% input - A is an N*N nonsingular matrix

- B is an N*1 matrix

- P is an N*1 matrix :the initial guess

- delta is the tolerance for P

- max1 is the maximum number of iterations

% output - X is an N*1 matrix : the jacobi approximation to the solution of AX=B

N=length(B);

for k=1:max1

for j=1:N

if j==1

X(1)=(B(1)-A(1,2:N)*P(2:N))/A(1,1);

elseif j==N

X(N)=(B(N)-A(N,1:N-1)*(X(1:N))')/A(N,N);

else

浙江工业大学化材学院 jackdong

X(j)=(B(j)-A(j,1:j-1)*X(1:j-1)'-A(j,j+1:N)*P(j+1:N)/A(j,j);

end

err=abs(norm(X'-P));

relerr=err/(norm(X)+eps);

P=X';

if (errbreak;

end end

X=X';

----------------------------------------------------------------------------------------------------------------------

5.3.非线性Seidel迭代

功能：求解非线性不动点方程组X=G(X),给定初始近似值P0,并生成序列{Pk},收敛到解P ----------------------------------------------- function [P,iter]=seidel(G,P,delta,max1)

% input - G is an N*N nonlinear matrix

% - P is the initial guess at solution

% - delta is the error bound

% - max1 is the maximum number of iterations

% output - P is the seidel approimation to the solution

% - iter is the number of iterations required

浙江工业大学化材学院 jackdong

N=length(P);

for k=1:max1

X=P;

for j=1:N

A=feval('G',X);

X(j)=A(j);

end

err=abs(norm(X-P));

relerr=err/(norm(X)+eps);

P=X;

iter=k;

if (errbreak;

end end

----------------------------------------------------------------------------------------------------------------------

5.4. 牛顿-拉夫森法

功能：求解非线性方程组F(X)=0，给定初始近似值P0，并生成序列{Pk},收敛到解P -----------------------------------------------

function [P,iter,err]=newdim(F,JF,P,delta,epsilon,max1)

% input - F is the system saved as the M-file F.m

% - JF is the Jacobian of F saved as the M-file JF.m

浙江工业大学化材学院 jackdong

% - P is the initial approximation to the solution

% - delta is the tolerance for P

% - max1 is the maximum number of iterations

% output - P is the approimation to the solution

% - iter is the number of iterations required

% - err is the error estimate for P

Y=feval(F,P);

for k=1:max1

J=feval(JF,P);

Q=P-(J\\Y);

Z=feval(F,Q);

err=norm(Q-P);

err=abs(norm(X-P));

relerr=err/(norm(Q)+eps);

P=Q;

Y=Z;

iter=k;

if (errbreak;

end end

----------------------------------------------------------------------------------------------------------------------

浙江工业大学化材学院 jackdong

5.5. 超松弛迭代

功能：用超松弛迭代法求解线性方程组

---------------------------------------------

function [x,k]=SORmethod(A,b,x0,N,emg,w)

% A:线性方程组左端矩阵

% b:线性方程组右端向量

% x0:迭代初值

% N:迭代次数上界，若迭代次数大于n，则迭代失败

% emg:精度指标

% w:松弛因子

% x:用迭代法求得的线性方程组的近似解

n=length(A);

x1=zeros(n,1);x2=zeros(n,1);

x1=x0;

r=max(abs(b-A*x1)); k=0;

while r>emg

for i=1:n

sum=0;

for j=1:n

if j>i

sum=sum+A(i,j)*x1(j);

浙江工业大学化材学院 jackdong

elseif jsum=sum+A(i,j)*x2(j);

end

x2(i)=x1(i)+w*(b(i)-sum)/A(i,i);

end

r=max(abs(x2-x1));

x1=x2;

k=k+1;

if k>N

disp('迭代失败，返回')；

return;

end end

x=x1;

----------------------------------------------------------------------------------------------------------------------

5.6. 对超松弛迭代

功能：用对超松弛(SSOR)迭代法求解线性方程组 --------------------------------------------------

function [x,k]=SSORmethod(A,b,x0,N,emg,w)

% A:线性方程组左端矩阵

% b:线性方程组右端向量

% x0:迭代初值

浙江工业大学化材学院 jackdong

% N:迭代次数上界，若迭代次数大于n，则迭代失败

% emg:精度指标

% w:松弛因子

% x:用迭代法求得的线性方程组的近似解

n=length(A);

x1=zeros(n,1);x2=zeros(n,1); x3=zeros(n,1);

x1=x0;

r=max(abs(b-A*x1)); k=0;

while r>emg

for i=1:n

sum=0;

for j=1:n

if j>i

sum=sum+A(i,j)*x1(j);

elseif jsum=sum+A(i,j)*x2(j);

end

x2(i)=(1-w)*x1(i)+w*(b(i)-sum)/A(i,i);

end

for i=n:-1:1

浙江工业大学化材学院 jackdong

sum=0;

for j=1:n

if j>i

sum=sum+A(i,j)*x3(j);

elseif jsum=sum+A(i,j)*x2(j);

end

x3(i)=(1-w)*x2(i)+w*(b(i)-sum)/A(i,i);

end

r=max(abs(x3-x1));

x1=x3;

k=k+1;

if k>N

disp('迭代失败，返回')；

return;

end end

x=x1;

浙江工业大学化材学院 jackdong

第六章线性方程组的直接法

6.1. 追赶法

功能：用追赶法解三对角线性方程组A*x=f,A为三对角矩阵

------------------------------------------------------- function x=threedia(a,b,c,f)

% a:三对角矩阵A的下对角线元素,a(1)=0

% b:三对角矩阵A的对角线元素

% c:三对角矩阵A的上对角线元素,c(N)=0

% f:方程组的右端向量

N=length(f);

x=zeros(1,N);y=zeros(1,N); d=zeros(1,N);u=zeros(1,N);

d(1)=b(1);

for i=1:N-1

u(i)=c(i)/d(i);

d(i+1)=b(i+1)-a(i+1)*u(i); end

% 追的过程

y(1)=f(1)/d(1);

for i=2:N

y(i)=(f(i)-a(i)*y(i-1))/d(i); end

浙江工业大学化材学院 jackdong

% 赶的过程

x(N)=y(N);

for i=N-1:-1:1

x(i)=y(i)-u(i)*x(i+1); end

----------------------------------------------------------------------------------------------------------------------

6.2. Cholesky方法

功能：用Cholesky分解法求解对称方程组A*x=b

----------------------------------------------------- function x=Chol_decompose(A,b)

% A:对称矩阵

% b:方程组的右端向量

% L:单位下三角矩阵

% D:单位上三角矩阵

% 对矩阵A进行三角分解：A=LDL'

N=length(A);

L=zeros(N,N);D=zeros(1,N);

for i=1:N

L(i,i)=1; end

D(1)=A(1,1);

for i=2:N

浙江工业大学化材学院 jackdong

for j=1:i-1

if j==1

L(i,j)=A(i,j)/D(j);

else

sum1=0;

for k=1:j-1

sum1=sum1+ L(i,k)*D(k)*L(j,k);

end

L(i,j)=(A(i,j)-sum1)/D(j);

end

sum2=0;

for k=1:i-1

sum2=sum2+L(i,k)^2*D(k);

end

D(i)=A(i,i)-sum2; end

% 分别求解线性方程组Ly=b;L'x=y/D

y=zeros(1,N);

y(1)=b(1);

for i=2:N

sumi=0;

浙江工业大学化材学院 jackdong

for k=1:i-1

sumi=sumi+L(i,k)*y(k);

end

y(i)=b(i)-sumi; end

x=zeros(1,N);

x(N)=y(N)/D(N);

for i=N-1:-1:1

sumi=0;

for k=i+1:N

sumi=sumi+L(k,i)*x(k);

end

x(i)=y(i)/D(i)-sumi; end

----------------------------------------------------------------------------------------------------------------------

6.3. 矩阵分解方法

功能：基于Gauss消去法的LU分解求解线性方程组A*x=b

------------------------------------------------------ function x=lu_decompose(A,b)

% A:系数矩阵

% b:方程组的右端向量

% L:单位下三角矩阵 n=length(b);

浙江工业大学化材学院 jackdong

L=eye(n);U=zeros(n,n);

x=zeros(n,1);y=zeros(n,1);

for i=1:n

U(1,i)=A(1,i);

if i==1

L(i,1)=1;

else

L(i,1)=A(i,1)/U(1,1);

end end

for i=2:n

for j=1:n

sum=0;

for k=1:i-1

sum=sum+ L(i,k)*U(k,j);

end

U(i,j)=A(i,j)-sum;

if j~=n

sum=0;

for k=1:i-1

sum=sum+L(j+1,k)*U(k,i);

end

浙江工业大学化材学院 jackdong

L(j+1,1)=(A(j+1,i)-sum)/U(i,i);

end

end end

% 求解线性方程组Ly=b

y(1)=b(1);

for i=2:n

sum=0;

for j=1:k-1

sum=sumi+L(k,j)*y(j);

end

y(k)=b(k)-sum; end

% 解方程组Ux=y

x(n)=y(n)/U(n,n);

for i=n-1:-1:1

sum=0;

for j=k+1:n

sum=sum+U(k,i)*x(j);

end

x(k)=(y(k)-sum)/U(k,k); end

----------------------------------------------------------------------------------------------------------------------

浙江工业大学化材学院 jackdong

6.4. 消去法

功能：用Gauss列主元素消去法求解线性方程组A*x=b

-------------------------------------------------------------- function x=Gauss_pivot(A,b)

% A:系数矩阵

% b:方程组的右端向量

n=length(b);

x=zeros(n,1);c=zeros(1,n);d1=0;

for i=1:n-1

max=abs(A(i,i));

m=i;

for j=i+1:n

if maxm=j;

end end

if m~=1

for k=i:n

c(k)=A(i,k);

A(i,k)=A(m,k);

A(m,k)=c(k);

end

d1=b(i);

浙江工业大学化材学院 jackdong

b(i)=b(m);

b(m)=d1;

end

for k=i+1:n

for j=i+1:n

A(k,j)=A(k,j)-A(i,j)*A(k,j)/A(i,i);

end

b(k)=b(k)-b(i)*A(k,j)/A(i,i);

A(k,i)=0;

end

x(n)=b(n)/A(n,n);

for i=n-1:-1:1

sum=0;

for j=i+1:n

sum=sum+ A(i,j)*x(j);

end

x(i)=(b(i)-sum)/A(i,i); end

======================================================================

浙江工业大学化材学院 jackdong

第七章数值优化

7.1. 黄金分割法求极小值

功能：用黄金分割法求f(x)在区间[a,b]上的近似极小值。当且仅当f(x)在[a,b]上为单峰时此方法适用

----------------------------------------------------------------- function [S,E,G]=golden(f,a,b,delta,epsilon)

% input - f is the object function input as a string 'f'

% - a and b are the end points of the interval

% - delta is the tolerance for the abscissas

% - epsilon is the tolerance for the ordinates

% output - S=(p,yp) contains the abscissa p and the ordinate yp of the minimum

% - E=(dp,dy) contains the error bounds for p and yp

% -G is an n*4 matrix : the kth row contains

% -[ak ck dk bk]: the values of a,c,d and b at the kth interation

r1= (sqrt(5)-1)/2;

r2=r1^2;

h=b-a;

ya=feval(f,a);

yb=feval(f,b);

c=a+r2*h;

d=a+r1*h;

yc=feval(f,c);

浙江工业大学化材学院 jackdong

yd=feval(f,d); k=1;

A(k)=a;B(k)=b;C(k)=c;D(k)=d;

while (abs(yb-ya)>epsilon) | (h>delta)

k=k+1;

if (ycb=d;

yb=yd;

d=c;

yd=yc;

h=b-a;

c=a+r2*h;

yc=feval(f,c);

else

a=c;

ya=yc;

c=d;

yc=yd;

h=b-a;

d=a+r1*h;

ya=feval(f,d);

end

浙江工业大学化材学院 jackdong

A(k)=a;B(k)=b;C(k)=c;D(k)=d; end

dp=abs(b-a);

dy=abs(yb-ya); p=a;

yp=ya;

if (ybp=b;

yp=yb; end;

G=[A' C' D' B'];

S=[p yp];

E=[dp dy];

----------------------------------------------------------------------------------------------------------------------

7.2. 斐波那契法求极小值

功能：用斐波那契法求f(x)在区间[a,b]上的近似极小值。当且仅当f(x)在[a,b]上为单峰时此方法适用

----------------------------------------------------------------- function X=fibonacci(f,a,b,tol,e)

% input - f is the object function input as a string 'f'

% - a and b are the end points of the interval

% - tol : length of uncertainy

% - e :distinguishability of constant

浙江工业大学化材学院 jackdong

% output - X : x and y coordinates of minimum

% note : this function calls the m-file fib.m while had been edited i=1; F=1;

while F<=(b-a)/tol

F=fib(i);

i=i+1; end

n=i-1;

A=zeros(1,n-2);B=zeros(1,n-2);

A(1)=a;

B(1)=b;

c=A(1)+(fib(n-2)/fib(n))*(B(1)-A(1));

d=A(1)+(fib(n-1)/fib(n))*(B(1)-A(1)); k=1;

while k=n-3

if feval(f,c)>feval(f,d)

A(k+1)=c;

B(k+1)=B(k);

c=d;

d=A(k+1)+(fib(n-k-1)/fib(n-k))*(B(k+1)-A(k+1));

浙江工业大学化材学院 jackdong

else

A(k+1)=A(k);

B(k+1)=d;

d=c;

c=A(k+1)+(fib(n-k-2)/fib(n-k))*(B(k+1)-A(k+1));

end

k=k+1; end

if feval(f,c)>feval(f,d)

A(n-2)=c;

B(n-2)=B(n-3);

c=d;

d=A(n-2)+(0.5+e)*(B(n-2)-A(n-2));

else

A(n-2)=A(N-3);

B(n-2)=d;

d=c;

c=A(n-2)+(0.5+e)*(B(n-2)-A(n-2)); end

if feval(f,c)>feval(f,d)

a=c;

b=B(n-2);

浙江工业大学化材学院 jackdong

else

a=A(n-2);

b=d; end

X=[(a+b)/2 feval(f,(a+b)/2)];

------------------------------------------------------------- % 下面的程序用来生成斐波那契数

function y=fib(n)

fz(1)=1;

fz(2)=1;

for k=3:n

fz(k)=fz(k-1)+fz(k-1); end

y=fz(n);

----------------------------------------------------------------------------------------------------------------------

7.3. 用2次插值求局部最小值

功能：求f(x)在区间[a,b]上的局部极小值，从初始近似点p0开始，然后在区间[a,p0]和[p0,b]上进行搜索

------------------------------------------------------------------- function [p,yp,dp,dy,P]=quadmin(f,a,b,delta,epsilon)

% input - f is the object function input as a string 'f'

% - a and b are the end points of the interval

% - delta is the tolerance for the abscissas

% - epsilon is the tolerance for the ordinates

浙江工业大学化材学院 jackdong

% output - p is the abscissa of the minimum

% - yp is the ordinate of the minimum

% - dp is the error bound for p

% - dy is the error bound for yp

% - P is the vector of interations

p0=a;

maxj=20;

maxk=30;

big=1e6;

err=1;

P(k)=p0;

cond=0; h=1;

if (abs(p0)>1e4),h=abs(p0)/1e4;end

while (kepsilon & cond~=5)

f1=feval(f,p0+0.00001)-feval(f,p0-0.00001)/0.00002;

if (f1>0),h=-abs(h);end

p1=p0+h;

p2=p0+2*h;

pmin=p0;

y0=feval(f,p0);

y1=feval(f,p1);

浙江工业大学化材学院 jackdong

y2=feval(f,p2);

ymin=y0;

cond=0;

j=0;

while (jdelta & cond==0)

if (y0<=y1)

p2=p1;

y2=y1;

h=h/2;

p1=p0+h;

y1=feval(f,p1);

else

if y2p1=p2;

y1=y2;

h=2*h;

p2=p0+2*h;

y2=feval(f,p2);

else

cond=-1;

end

浙江工业大学化材学院 jackdong

j=j+1;

if (abs(h)>big | abs(p0)>big), cond=5;end

end

if cond==5

pmin=p1;

ymin=feval(f,p1);

else

d=4*y1-2*y0-2*y2;

if d<0

hmin=h*(4*y1-3*y0-y2)/d;

else

hmin=h/3;

cond=4; end

pmin=p0+hmin;

ymin=feval(f,pmin);

h=abs(h);

h0=abs(hmin);

h1=abs(hmin-h);

h2=abs(hmin-2*h);

if (h0if (h179

浙江工业大学化材学院 jackdong

if (h2if (h==0),h=hmin;end

if (hif (abs(h)>big) | abs(pmin>big),cond=5;end

e0=abs(y0-ymin);

e1=abs(y1-ymin);

e2=abs(y2-ymin);

if (e0~=0 & e0if (e1~=0 & e1if (e2~=0 & e2if (e0~=0 & e1==0 & e2==0),err==0;end

if (errp0=pmin;

k=k+1;

P(k)=P0; end

if (cond==2 & hp=p0;

dp=h;

yp=feval(f,p);

dy=err;

----------------------------------------------------------------------------------------------------------------------

浙江工业大学化材学院 jackdong

7.4. 内德-米德法求最小值

功能：求f(x1,x2,......xn)的近似最小值，其中f是N个实量的连续函数，给定N+1个初始点Vk=(vk1,vk2......vkN),k=0,1,......N

---------------------------------------------------

function [V0,y0,dV,dy]=nelder(F,V,min1,max1,epsilon,show)

% input - F is the object function input as a string 'F'

- V is a 3*n matrix containing staring simplex

- min1 & max1 are minimum and maximum number of interations

- epsilon is the tolerance

- show==1 displays interations (P and Q)

% output -V0 is the vertex for the minimum

- y0 is the function value F(V0)

- dV is the size of the final simplex

- dy is the error bound for the minimum

- P is a matrix containing the vertex iteration

- Q is an array containing the interations for F(P)

if nargin ==5

show=0; end

[mm n]=size(V);

for j=1:n+1

Z=V(j,1:n);

Y(j)=feval(F,Z);

浙江工业大学化材学院 jackdong

end

[mm lo]=min(Y);

[mm hi]=max(Y);

li=hi;

ho=lo;

for j=1:n+1

if (j~=lo & j~=ho & j<=Y(li))

li=j;

end

if (j~=hi & j~=lo & Y(j)>=Y(ho))

ho=j;

end end

cnt=0;

while (Y(hi)>Y(lo)+epsilon & cntS=zeros(1,n+1);

for j=1:n+1

S=S+V(j,1:n);

end

M=(S-V(hi,1:n))/n;

R=2*M-V(hi,1:n);

yR=feval(F,R);

浙江工业大学化材学院 jackdong

if (yRif Y((li)V(hi,1:n)=R;

Y(hi)=yR;

else

E=2*R-M;

yE=feval(F,E);

if (yEV(hi,1:n)=E;

Y(hi)=yE;

else

V(hi,1:n)=R;

Y(hi)=yR;

end

else

if (yRV(hi,1:n)=R;

Y(hi)=yR;

end

C=(V(hi,1:n)+M)/2;

yC=feval(F,C);

浙江工业大学化材学院 jackdong

C2=(M+R)/2;

yC2=feval(F,C2);

if (yC2C=C2;

yC=yC2;

end

if (yCV(hi,1:n)=C;

Y(hi)=yC;

else

for j=1:n+1

if (j~=lo)

V(j,1:n)=(V(j,1:n)+V(lo,1:n))/2;

Z=V(j,1:n);

Y(j)=feval(F,Z);

end

[mm lo]=min(Y);

[mm hi]=max(Y);

li=hi;

ho=lo;

浙江工业大学化材学院 jackdong

for j=1:n+1

if (j~=lo & j~=hi & Y(j)<=Y(li))

li=j;

end

if (j~=hi & j~=lo & Y(j)>=Y(ho))

ho=j;

end end

cnt=cnt+1;

P(cnt,:)=V(lo,:);

Q(cnt)=Y(lo); end

snorm=0;

for j=1:n+1

s=norm(V(j)-V(lo));

if (s>=snorm)

snorm=s;

end end

Q=Q';

V0=V(lo,1:n)

y0=Y(lo);

浙江工业大学化材学院 jackdong

dV=snorm;

dy=abs(Y(hi)-Y(lo));

if (show==1)

disp(P);

disp(Q); end

----------------------------------------------------------------------------------------------------------------------

7.5. 最速下降法或梯度法

功能：用数值方法求f(X)的近似局部极小值，其中f为N个实数变量的连续函数，X=(x1,x2,,.....xN),从P0开始，适用梯度方法

----------------------------------------------------------

function [P0,y0,err]=grads(F,G,P0,max1,delta,epsilon,show)

% input - F is the object function input as a string 'F'

- G=-(i/norm(grad F))*grad F: the search direction input as a string 'G'

- P0 is the initial starting point

- max1 is the maximum number of interations

- delta is the tolerance for hmin in the single parameter minimization in the search direction

- epsilon is the tolerance for the err in y0

- show : if show==1 the interations are displayed

% output - P0 is the point for the minimum

- y0 is the function value F(P0)

- err is the error bound for y0

浙江工业大学化材学院 jackdong

- P is a vector containing the iterations

if nargin==5,show=0;end

[mm n]=size(P0);

maxj=10;

big=1e8 h=1;

P=zeros(maxj,n+1);

len=norm(P0);

y0=feval(F,P0);

if (len>1e4),h=len/1e4;end

cond=0;err=1;cnt=0;

P(cnt+1,:)=[P0 y0];

while (cntdelta | err>epsilon))

S=feval(G,P0);

P1=P0+h*S;

P2=P0+2*h*S;

y1=feval(F,P1);

y2=feval(F,P2);

cond=0;

j=0;

while (jlen=norm(P0)

浙江工业大学化材学院 jackdong

if (y0P2=P1;

y2=y1;

h=h/2;

P1=P0+h*S;

y1=feval(F,P1);

else

if y2P1=P2;

y1=y2;

h=2*h;

P2=P0+2*h*S;

y2=feval(F,P2);

else

cond=-1;

end

j=j+1;

if (hbig | len>big),cond=5;end

end

if cond==5

浙江工业大学化材学院 jackdong

Pmin=P1;

ymin=y1;

else

d=4*y1-2*y0-2*y2;

if d<0

hmin=h*(4*y1-3*y0-y2)/d;

else

cond=4;

hmin=h/3;

end

Pmin=P0+hmin*S;

ymin=feval(F,Pmin);

h=abs(h);

h0=abs(hmin);

h1=abs(hmin-h);

h2=abs(hmin-2*h);

if (h0if (h1if (h2if (h==0),h=hmin;end

if (he0=abs(y0-ymin);

浙江工业大学化材学院 jackdong

e1=abs(y1-ymin);

e2=abs(y2-ymin);

if (e0~=0 & e0if (e1~=0 & e1if (e2~=0 & e2if (e0==0 & e1==0 & e2==0),err==0;end

if (errif (cond==2 & hcnt=cnt+1;

P(cnt+1,:)=[Pmin ymin];

P0=Pmin;

y0=ymin; end

if show==1

disp(P)； end

======================================================================

因篇幅问题不能全部显示，请点此查看更多更全内容

查看全文

全部栏目

Matlab数值计算程序