数值计算方法matlab程序

function [x0,k]=bisect1(fun1,a,b,ep)

if nargin<4

ep=1e-5;

end

fa=feval(fun1,a);

fb=feval(fun1,b);

if fa*fb>0

x0=[fa,fb];

k=0;

return;

end

k=1;

while abs(b-a)/2>ep

x=(a+b)/2;

fx=feval(fun1,x);

if fx*fa<0

b=x;

fb=fx;

else

a=x;

fa=fx;

k=k+1;

end

end

x0=(a+b)/2;

>> fun1=inline('x^3-x-1');

>> [x0,k]=bisect1(fun1,1.3,1.4,1e-4)

x0 =

1.3247

k =

7

>>

简单迭代法

function [x0,k]=iterate1(fun1,x0,ep,N)

if nargin<4

N=500;

end

if nargin<3

ep=1e-5;

end

x=x0;

x0=x+2*ep;

while abs(x-x0)>ep & k

x0=x;

x=feval(fun1,x0);

k=k+1;

end

x0=x;

if k==N

warning('已达最大迭代次数')

end

>> fun1=inline('(x+1)^(1/3)');

>> [x0,k]=iterate1(fun1,1.5)

x0 =

1.3247

k =

7

>> fun1=inline('x^3-1');

>> [x0,k]=iterate1(fun1,1.5)

x0 =

Inf

k =

9

>>

Steffesen加速迭代（简单迭代法的加速）function [x0,k]=steffesen1(fun1,x0,ep,N)

if nargin<4

N=500;

end

if nargin<3

ep=1e-5;

end

x=x0;

x0=x+2*ep;

k=0;

while abs(x-x0)>ep & k

x0=x;

y=feval(fun1,x0);

z=feval(fun1,y);

x=x0-(y-x0)^2/(z-2*y+x0);

k=k+1;

end

x0=x;

if k==N

warning('已达最大迭代次数')

end

>> fun1=inline('(x+1)^(1/3)');

>> [x0,k]=steffesen1(fun1,1.5)

x0 =

1.3247

k =

3

>> fun1=inline('x^3-1');

>> [x0,k]=steffesen1(fun1,1.5)

x0 =

1.3247

k =

6

Newton迭代

function [x0,k]=Newton7(fname,dfname,x0,ep,N)

if nargin<5

N=500;

end

if nargin<4

ep=1e-5;

end

x=x0;

x0=x+2*ep;

k=0;

while abs(x-x0)>ep & k

x0=x;

x=x0-feval(fname,x0)/feval(dfname,x0);

k=k+1;

end

x0=x;

if k==N

warning('已达最大迭代次数')

end

>> fname=inline('x-cos(x)');

>> dfname=inline('1+sin(x)');

>> [x0,k]=Newton7(fname,dfname,pi/4,1e-8)

x0 =

0.7391

k =4

非线性方程求根的Matlab函数调用举例：

1.求多项式的根：求f(x)=x^3-x-1=0的根：

>> roots([1 0 -1 -1])

ans =

1.3247

-0.6624 + 0.5623i

-0.6624 - 0.5623i

2.求一般函数的根

>> fun=inline('x*sin(x^2-x-1)','x')

fun =

Inline function:

fun(x) = x*sin(x^2-x-1)

>> fplot(fun,[-2 0.1]);grid on

>> x=fzero(fun,[-2,-1])

x =

-1.5956

>> x=fzero(fun,[-1 -0.1])

x =

-0.6180

[x,f,h]=fsolve(fun,-1.6)

x =

-1.5956

f =

1.4909e-009

h =

1

（h>0表示收敛，h<0表示发散，h=0表示已达到设定的计算函数值的最大次数）

第三章：线性方程组的数值解法

1. 高斯消元法

function [A,x]=gauss3(A,b)

%本算法用顺序高斯消元法求解线性方程组

n=length(b);

A=[A,b];

for k=1:n-1

A((k+1):n,(k+1):(n+1))=A((k+1):n,(k+1):(n+1))-A((k+1):n,k)/A(k,k)*A(k,(k+1):(n+1));

A((k+1):n,k)=zeros(n-k,1);

A;

end

x=zeros(n,1);

%上面为消元过程

x(n)=A(n,n+1)/A(n,n);

for k=n-1:-1:1

x(k)=(A(k,n+1)-A(k,(k+1):n)*x((k+1:n)))/A(k,k);

end

%上面为回代过程

>> A=[2 3 4;3 5 2;4 3 30];

>> b=[6,5,32]'

b =

6

5

32

>> [A,x]=gauss3(A,b)

A =

2.0000

3.0000

4.0000 6.0000

0 0.5000 -4.0000 -4.0000

0 0 -2.0000 -4.0000

x =

-13

8

2

列选主元的高斯消元法：

function [A,x]=gauss5(A,b)

%本算法用列选主元的高斯消元法求解线性方程组

n=length(b);

A=[A,b];

for k=1:n-1

%选主元

[ap,p]=max(abs(A(k:n,k)));

p=p+k-1;

if p>k

t=A(k,:);

A(k,:)=A(p,:);

A(p,:)=t;

end

%消元

A((k+1):n,(k+1):(n+1))=A((k+1):n,(k+1):(n+1))-A((k+1):n,k)/A(k,k)*A(k,(k+1):(n+1));

A((k+1):n,k)=zeros(n-k,1);

end

%回代

x=zeros(n,1);

x(n)=A(n,n+1)/A(n,n);

for k=n-1:-1:1

x(k)=(A(k,n+1)-A(k,(k+1):n)*x((sk+1:n)))/A(k,k);

end

>> A=[2 3 4;3 5 2;4 3 30]; b=[6,5,32]';

>> [A,x]=gauss5(A,b)

A =

4.0000 3.0000 30.0000 32.0000

0 2.7500 -20.5000 -19.0000

0 0 0.1818 0.3636

x =

-13

8

2

三角分解法：Doolittle 分解

function [L,U]=doolittle1(A)

n=length(A);

U=zeros(n);

L=eye(n);

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

L(2:n,1)=A(2:n,1)/U(1,1);

for k=2:n

U(k,k:n)=A(k,k:n)-L(k,1:k-1)*U(1:k-1,k:n);

L(k+1:n,k)=A(k+1:n,k)-L(k+1:n,1:k-1)*U(1:k-1,n)/U(k,k); End

y=zeros(n,1);

x=y;

y(1)=b(1);

for i=2:n

y(i)=b(i)-L(i,1:i-1)*y(1:i-1);

end

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

for i=n-1:-1:1

x(i)=(y(i)-U(i,i+1:n)*x(i+1:n))/U(i,i);

end

>> A=[1 2 3;2 5 2 ;3 1 5];b=[14 18 20]';

>> [L,U,x]=doolittle1(A,b)

L =

1 0 0

2 1 0

3 -8 1

U =

1 2 3

0 1 -4

0 0 -36

x =

2.8333

1.3333

2.8333

平方根法：

function [L,x]=choesky3(A,b)

n=length(A);

L=zeros(n);

L(:,1)=A(:,1)/sqrt(A(1,1));

for k=2:n

L(k,k)=A(k,k)-L(k,1:k-1)*L(k,1:k-1)';

L(k,k)=sqrt(L(k,k));

for i=k+1:n

L(i,k)=(A(i,k)-L(i,1:k-1)*L(k,1:k-1)')/L(k,k);

end

end

y=zeros(n,1);

x=y;

y(1)=b(1)/L(1,1);

for i=2:n

y(i)=(b(i)-L(i,1:i-1)*y(1:i-1))/L(i,i);

end

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

for i=n-1:-1:1

x(i)=(y(i)-L(i+1:n,i)'*x(i+1:n))/L(i,i);

end

>> A=[4 -1 1;-1 4.25 2.75;1 2.75 3.5]

A =

4.0000 -1.0000 1.0000

-1.0000 4.2500 2.7500

1.0000

2.7500

3.5000

>> b=[4 6 7.25]'

b =

4.0000

6.0000

7.2500

[L,x]=choesky3(A,b)

L =

2.0000 0 0

-0.5000 2.0000 0

0.5000 1.5000 1.0000

x =

1

1

1

>>

迭代法求方程组的解

Jacobi迭代法：

function [x,k]=jacobi2(a,b,x0,ep,N)

%本算法用Jacobi迭代求解ax=b,用分量形式n=length(b);

k=0;

if nargin<5

N=500;

end

if nargin<4

ep=1e-5;

end

if nargin<3

x0=zeros(n,1);

y=zeros(n,1);

end

x=x0;x0=x+2*ep;

while norm(x-x0,inf)>ep & k

k=k+1;

x0=x;

for i=1:n

y(i)=b(i);

for j=1:n

if j~=i

y(i)=y(i)-a(i,j)*x0(j);

end

end

if abs(a(i,i))<1e-10|k==N

warning('a(i,i) is too small');

return

end

y(i)=y(i)/a(i,i);

end

x=y;

end

a=[4 3 0;3 4 -1; 0 -1 4];b=[24 30 -24]'；

[x,k]=jacobi2(a,b)

x =

3.0000

4.0000

-5.0000

k =

59

Gauss-seidel迭代法：

function [x,k]=gaussseide2(a,b,x0,ep,N)

%本算法用Gauss-seidel迭代求解ax=b,用分量形式n=length(b);

k=0;

if nargin<5

N=500;

end

if nargin<4

ep=1e-5;

end

if nargin<3

x0=zeros(n,1);

y=zeros(n,1);

end

x=x0;x0=x+2*ep;

while norm(x-x0,inf)>ep & k

k=k+1;

x0=x;

y=x;

for i=1:n

z(i)=b(i);

for j=1:n

if j~=i

z(i)=z(i)-a(i,j)*x(j);

end

end

if abs(a(i,i))<1e-10|k==N

warning('a(i,i) is too small');

return

end

z(i)=z(i)/a(i,i);

x(i)=z(i);

end

end

[x,k]=gaussseide2(a,b)

x =

3.0000

4.0000

-5.0000

k =

25

最速下降法function [x,k]=zuisuxiajiang(A,b,x0,ep,N)

%本算法用最速下降算法求解正定方程组Ax=b,

n=length(b);

if nargin<5

N=500;

end

if nargin<4

ep=1e-8;

end

if nargin<3

x0=ones(n,1);

end

x=x0;

x0=x+2*ep;

r=b-A*x;

d=r;

k=0;

while norm(x-x0,inf)>ep & k

k=k+1;

x0=x;

lamda=(d'*d)/(d'*A*d);

x=x0+lamda*d;

r=b-A*x;

d=r;

end

if k==N

warning('已达最大迭代次数')

end

共轭梯度算法function [x,k]=gongertidufa(A,b,x0,ep,N)

%本算法用共轭梯度算法求解正定方程组Ax=b,,

n=length(b);

if nargin<5

N=500;

end

if nargin<4

ep=1e-8;

end

if nargin<3

x0=ones(n,1);

end

x=x0;

x0=x+2*ep;

r=b-A*x;

d=r;

k=0;

while norm(x-x0,inf)>ep & k

x0=x;

lamda=(r'*r)/(d'*A*d);

r1=r;

x=x0+lamda*d;

r=b-A*x;

beta=(r'*r)/(r1'*r1);

d=r+beta*d;

end

if k==N

warning('已达最大迭代次数') end

常微分方程数值解

function [x,y]=Euler1(fun,xspan,y0,h)

%本算法用欧拉格式计算微分方程y'=f(x,y)的解。

%fun为右端函数，xspan为求解的自变量区间，yo为初值，h为步长。% 输出向量x以及对应的函数值向量y.

x=xspan(1):h:xspan(2);

y(1)=y0;

n=length(x);

for k=1:n-1

y(k+1)=y(k)+h*feval(fun,x(k),y(k));

end

x=x';y=y';

fun=inline('y-2*x/y');

[x,y1]=Euler1(fun,[0,1],1,0.1)

x =

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

y1 =

1

1.1

1.1918

1.2774

1.3582

1.4351

1.509

1.5803

1.6498

1.7178

1.7848

function [x,y]=Euler2(fun,xspan,y0,h)

%本算法用改进的欧拉格式计算微分方程y'=f(x,y)的解。

%fun为右端函数，xspan为求解的自变量区间，yo为初值，h为步长。% 输出向量x以及对应的函数值向量y.

x=xspan(1):h:xspan(2);

y(1)=y0;

n=length(x);

for k=1:n-1

k1=feval(fun,x(k),y(k));

y(k+1)=y(k)+h*k1;

k2=feval(fun,x(k+1),y(k+1));

y(k+1)=y(k)+h*(k1+k2)/2;

end

x=x';y=y';

[x,y2]=Euler2(fun,[0,1],1,0.1)

x =

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

y2 =

1

1.0959

1.1841

1.2662

1.3434

1.4164

1.486

1.5525

1.6165

1.6782

1.7379

y3=sqrt(1+2*x)

y3 =

1

1.0954

1.1832

1.2649

1.3416

1.4142

1.4832

1.5492

1.6125

1.6733

1.7321

plot(x,y1,'+',x,y2,'o',x,y3,'r')

四阶Runge-Kutta方法求解初值问题

function [x,y]=rungekutta(dfun,xspan,y0,h)

%本算法用四阶runge-kutta 方法求解初值问题y’=f(x,y),y(x0)=y0 x=xspan(1):h:xspan(2);

y(1)=y0;

N=length(x);

for n=1:N-1

k1=feval(dfun,x(n),y(n));

k2=feval(dfun,x(n)+h/2,y(n)+h/2*k1);

k3=feval(dfun,x(n)+h/2,y(n)+h/2*k2);

k4=feval(dfun,x(n+1),y(n)+h*k3);

y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6;

end

x=x';y=y';

>> dfun=inline('y-2*x/y');

>> [x,y]=rungekutta(dfun,[0,1],1,0.2)

x =

0.2

0.4

0.6

0.8

1

y =

1

1.1832

1.3417

1.4833

1.6125

1.7321

>> y-y3

ans =

1.3331e-005

2.6143e-005

4.1761e-005

6.2492e-005

9.1075e-005

>>