数值计算方法matlab程序.doc
《数值计算方法matlab程序.doc》由会员分享,可在线阅读,更多相关《数值计算方法matlab程序.doc(16页珍藏版)》请在冰点文库上搜索。
数值计算方法matlab程序
二分法
function[x0,k]=bisect1(fun1,a,b,ep)
ifnargin<4
ep=1e-5;
end
fa=feval(fun1,a);
fb=feval(fun1,b);
iffa*fb>0
x0=[fa,fb];
k=0;
return;
end
k=1;
whileabs(b-a)/2>ep
x=(a+b)/2;
fx=feval(fun1,x);
iffx*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)
ifnargin<4
N=500;
end
ifnargin<3
ep=1e-5;
end
x=x0;
x0=x+2*ep;
k=0;
whileabs(x-x0)>ep&kx0=x;
x=feval(fun1,x0);
k=k+1;
end
x0=x;
ifk==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)
ifnargin<4
N=500;
end
ifnargin<3
ep=1e-5;
end
x=x0;
x0=x+2*ep;
k=0;
whileabs(x-x0)>ep&kx0=x;
y=feval(fun1,x0);
z=feval(fun1,y);
x=x0-(y-x0)^2/(z-2*y+x0);
k=k+1;
end
x0=x;
ifk==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)
ifnargin<5
N=500;
end
ifnargin<4
ep=1e-5;
end
x=x0;
x0=x+2*ep;
k=0;
whileabs(x-x0)>ep&kx0=x;
x=x0-feval(fname,x0)/feval(dfname,x0);
k=k+1;
end
x0=x;
ifk==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([10-1-1])
ans=
1.3247
-0.6624+0.5623i
-0.6624-0.5623i
2.求一般函数的根
>>fun=inline('x*sin(x^2-x-1)','x')
fun=
Inlinefunction:
fun(x)=x*sin(x^2-x-1)
>>fplot(fun,[-20.1]);gridon
>>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];
fork=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);
fork=n-1:
-1:
1
x(k)=(A(k,n+1)-A(k,(k+1):
n)*x((k+1:
n)))/A(k,k);
end
%上面为回代过程
>>A=[234;352;4330];
>>b=[6,5,32]'
b=
6
5
32
>>[A,x]=gauss3(A,b)
A=
2.00003.00004.00006.0000
00.5000-4.0000-4.0000
00-2.0000-4.0000
x=
-13
8
2
列选主元的高斯消元法:
function[A,x]=gauss5(A,b)
%本算法用列选主元的高斯消元法求解线性方程组
n=length(b);
A=[A,b];
fork=1:
n-1
%选主元
[ap,p]=max(abs(A(k:
n,k)));
p=p+k-1;
ifp>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);
fork=n-1:
-1:
1
x(k)=(A(k,n+1)-A(k,(k+1):
n)*x((sk+1:
n)))/A(k,k);
end
>>A=[234;352;4330];b=[6,5,32]';
>>[A,x]=gauss5(A,b)
A=
4.00003.000030.000032.0000
02.7500-20.5000-19.0000
000.18180.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);
fork=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);
fori=2:
n
y(i)=b(i)-L(i,1:
i-1)*y(1:
i-1);
end
x(n)=y(n)/U(n,n);
fori=n-1:
-1:
1
x(i)=(y(i)-U(i,i+1:
n)*x(i+1:
n))/U(i,i);
end
>>A=[123;252;315];b=[141820]';
>>[L,U,x]=doolittle1(A,b)
L=
100
210
3-81
U=
123
01-4
00-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));
fork=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));
fori=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);
fori=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);
fori=n-1:
-1:
1
x(i)=(y(i)-L(i+1:
n,i)'*x(i+1:
n))/L(i,i);
end
>>A=[4-11;-14.252.75;12.753.5]
A=
4.0000-1.00001.0000
-1.00004.25002.7500
1.00002.75003.5000
>>b=[467.25]'
b=
4.0000
6.0000
7.2500
[L,x]=choesky3(A,b)
L=
2.000000
-0.50002.00000
0.50001.50001.0000
x=
1
1
1
>>
迭代法求方程组的解
Jacobi迭代法:
function[x,k]=jacobi2(a,b,x0,ep,N)
%本算法用Jacobi迭代求解ax=b,用分量形式
n=length(b);
k=0;
ifnargin<5
N=500;
end
ifnargin<4
ep=1e-5;
end
ifnargin<3
x0=zeros(n,1);
y=zeros(n,1);
end
x=x0;x0=x+2*ep;
whilenorm(x-x0,inf)>ep&kk=k+1;
x0=x;
fori=1:
n
y(i)=b(i);
forj=1:
n
ifj~=i
y(i)=y(i)-a(i,j)*x0(j);
end
end
ifabs(a(i,i))<1e-10|k==N
warning('a(i,i)istoosmall');
return
end
y(i)=y(i)/a(i,i);
end
x=y;
end
a=[430;34-1;0-14];b=[2430-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;
ifnargin<5
N=500;
end
ifnargin<4
ep=1e-5;
end
ifnargin<3
x0=zeros(n,1);
y=zeros(n,1);
end
x=x0;x0=x+2*ep;
whilenorm(x-x0,inf)>ep&kk=k+1;
x0=x;
y=x;
fori=1:
n
z(i)=b(i);
forj=1:
n
ifj~=i
z(i)=z(i)-a(i,j)*x(j);
end
end
ifabs(a(i,i))<1e-10|k==N
warning('a(i,i)istoosmall');
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);
ifnargin<5
N=500;
end
ifnargin<4
ep=1e-8;
end
ifnargin<3
x0=ones(n,1);
end
x=x0;
x0=x+2*ep;
r=b-A*x;
d=r;
k=0;
whilenorm(x-x0,inf)>ep&kk=k+1;
x0=x;
lamda=(d'*d)/(d'*A*d);
x=x0+lamda*d;
r=b-A*x;
d=r;
end
ifk==N
warning('已达最大迭代次数')
end
共轭梯度算法
function[x,k]=gongertidufa(A,b,x0,ep,N)
%本算法用共轭梯度算法求解正定方程组Ax=b,,
n=length(b);
ifnargin<5
N=500;
end
ifnargin<4
ep=1e-8;
end
ifnargin<3
x0=ones(n,1);
end
x=x0;
x0=x+2*ep;
r=b-A*x;
d=r;
k=0;
whilenorm(x-x0,inf)>ep&kk=k+1;
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
ifk==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);
fork=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
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);
fork=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
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);
forn=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
0.2
0.4
0.6
0.8
1
y=
1
1.1832
1.3417
1.4833
1.6125
1.7321
>>y-y3
ans=
0
1.3331e-005
2.6143e-005
4.1761e-005
6.2492e-005
9.1075e-005
>>