实验一非线性方程的数值解法(一)信息与计算科学金融崔振威201002034031
一、实验目的:
熟悉二分法和简单迭代法的算法实现。
二、实验内容:
教材P40 2.1.5
三、实验要求
1 根据实验内容编写二分法和简单迭代法的算法实现
2 简单比较分析两种算法的误差
3 试构造不同的迭代格式,分析比较其收敛性
(一)、二分法程序:
function ef=bisect(fx,xa,xb,n,delta)
% fx是由方程转化的关于x的函数,有fx=0。
% xa 解区间上限
% xb 解区间下限
% n 最多循环步数,防止死循环。
%delta 为允许误差
x=xa;fa=eval(fx);
x=xb;fb=eval(fx);
disp(' [ n xa xb xc fc ]');
for i=1:n
xc=(xa+xb)/2;x=xc;fc=eval(fx);
X=[i,xa,xb,xc,fc];
disp(X),
if fc*fa<0
xb=xc;
else xa=xc;
end
if (xb-xa) end (二)、简单迭代法程序: function [x0,k]=iterate (f,x0,eps,N) if nargin<4 N=500; end if nargin<3 ep=1e-12; end x=x0; x0=x+2*eps; k=0; while abs(x-x0)>eps & k x0=x; x=feval(f,x0); k=k+1; end x0=x; if k==N end 解:a、g(x)=x5-3x3-2x2+2 二分法求方程: (1)、在matlab的命令窗口中输入命令: >> fplot('[x^5-3*x^3-2*x^2+2]',[-3,3]);grid 得下图: 由上图可得知:方程在[-3,3]区间有根。 (2)、二分法输出结果 >> f='x^5-3*x^3-2*x^2+2' f = x^5-3*x^3-2*x^2+2 >> bisect(f,-3,3,20,10^(-12)) 2.0000 - 3.0000 0 -1.5000 0.0313 3.0000 -3.0000 -1.5000 -2.2500 -31.6182 4.0000 -2.2500 -1.5000 -1.8750 -8.4301 5.0000 -1.8750 -1.5000 -1.6875 -2.9632 6.0000 -1.6875 -1.5000 -1.5938 -1.2181 7.0000 -1.5938 -1.5000 -1.5469 -0.5382 8.0000 -1.5469 -1.5000 -1.5234 -0.2405 9.0000 -1.5234 -1.5000 -1.5117 -0.1015 10.0000 -1.5117 -1.5000 -1.5059 -0.0343 11.0000 -1.5059 -1.5000 -1.5029 -0.0014 12.0000 -1.5029 -1.5000 -1.5015 0.0150 13.0000 -1.5029 -1.5015 -1.5022 0.0068 14.0000 -1.5029 -1.5022 -1.5026 0.0027 15.0000 -1.5029 -1.5026 -1.5027 0.0007 16.0000 -1.5029 -1.5027 -1.5028 -0.0003 17.0000 -1.5028 -1.5027 -1.5028 0.0002 18.0000 -1.5028 -1.5028 -1.5028 -0.0001 19.0000 -1.5028 -1.5028 -1.5028 0.0001 20.0000 -1.5028 -1.5028 -1.5028 -0.0000 2、迭代法求方程: 迭代法输出结果: >> f=inline('x^5-3*x^3-2*x^2+2'); >> [x0,k]=iterate(fun1,2) x0 = 2 k = 1 >> [x0,k]=iterate(fun1,1.5) x0 = NaN k = 6 >> [x0,k]=iterate(fun1,2.5) x0 = NaN k = 5 (3)、误差分析:由二分法和迭代法输出结果可知,通过定点迭代法得出方程的解误差比二分法大,而利用二分法求出的结果中,可以清楚看出方程等于零时的解,其误差比迭代法小。 b、g(x)=cos(sin(x)) 二分法求方程: (1)、在matlab的命令窗口中输入命令: >> fplot('[cos(sin(x))]',[-4,4]);grid 得下图: 由上图可得知:方程在[-4,4]区间无根。 (2)、二分法输出结果 >>f='cos(sin(x))' f = cos(sin(x)) >> bisect(f,-4,4,20,10^(-12)) 2.0000 0 4.0000 2.0000 0.6143 3.0000 2.0000 4.0000 3.0000 0.9901 4.0000 3.0000 4.0000 3.5000 0.9391 5.0000 3.5000 4.0000 3.7500 0.8411 6.0000 3.7500 4.0000 3.8750 0.7842 7.0000 3.8750 4.0000 3.9375 0.7554 8.0000 3.9375 4.0000 3.9688 0.7412 9.0000 3.9688 4.0000 3.9844 0.7341 10.0000 3.9844 4.0000 3.9922 0.7305 11.0000 3.9922 4.0000 3.9961 0.7288 12.0000 3.9961 4.0000 3.9980 0.7279 13.0000 3.9980 4.0000 3.9990 0.7275 14.0000 3.9990 4.0000 3.9995 0.7273 15.0000 3.9995 4.0000 3.9998 0.7271 16.0000 3.9998 4.0000 3.9999 0.7271 17.0000 3.9999 4.0000 3.9999 0.7271 18.0000 3.9999 4.0000 4.0000 0.7270 19.0000 4.0000 4.0000 4.0000 0.7270 20.0000 4.0000 4.0000 4.0000 0.7270 2、迭代法求方程: 迭代法输出结果: >> f=inline('cos(sin(x))'); >> [x0,k]=iterate(f,0.5) x0 = 0.7682 k = 15 >> [x0,k]=iterate(f,1) x0 = 0.7682 k = 15 >> [x0,k]=iterate(f,1.5) x0 = 0.7682 k = 16 >> [x0,k]=iterate(f,2) x0 = 0.7682 k = 15 >> [x0,k]=iterate(f,2.5) x0 = 0.7682 k = 14 (3)、由于该方程无解,所以无法比较误差。 c、g(x)=x2-sin(x+0.15) 二分法求方程: (1)、在matlab的命令窗口中输入命令: >> fplot('[x^2-sin(x+0.15)]',[-10,10]);grid 得下图: 由上图可得知:方程在[-3,3]区间有根。 (2)、二分法输出结果 >> f='x^2-sin(x+0.15)' f = x^2-sin(x+0.15) >> bisect(f,-3,3,30,10^(-12)) 1.0000 -3.0000 3.0000 0 -0.1494 2.0000 - 3.0000 0 -1.5000 3.2257 3.0000 -1.5000 0 -0.7500 1.1271 4.0000 -0.7500 0 -0.3750 0.3637 5.0000 -0.3750 0 -0.1875 0.0726 6.0000 -0.1875 0 -0.0938 -0.0474 7.0000 -0.1875 -0.0938 -0.1406 0.0104 8.0000 -0.1406 -0.0938 -0.1172 -0.0191 9.0000 -0.1406 -0.1172 -0.1289 -0.0045 10.0000 -0.1406 -0.1289 -0.1348 0.0029 11.0000 -0.1348 -0.1289 -0.1318 -0.0008 12.0000 -0.1348 -0.1318 -0.1333 0.0011 13.0000 -0.1333 -0.1318 -0.1326 0.0001 14.0000 -0.1326 -0.1318 -0.1322 -0.0003 15.0000 -0.1326 -0.1322 -0.1324 -0.0001 16.0000 -0.1326 -0.1324 -0.1325 0.0000 17.0000 -0.1325 -0.1324 -0.1324 -0.0000 18.0000 -0.1325 -0.1324 -0.1325 -0.0000 19.0000 -0.1325 -0.1325 -0.1325 0.0000 20.0000 -0.1325 -0.1325 -0.1325 0.0000 21.0000 -0.1325 -0.1325 -0.1325 0.0000 22.0000 -0.1325 -0.1325 -0.1325 0.0000 23.0000 -0.1325 -0.1325 -0.1325 -0.0000 24.0000 -0.1325 -0.1325 -0.1325 0.0000 25.0000 -0.1325 -0.1325 -0.1325 -0.0000 26.0000 -0.1325 -0.1325 -0.1325 0.0000 27.0000 -0.1325 -0.1325 -0.1325 0.0000 28.0000 -0.1325 -0.1325 -0.1325 0.0000 29.0000 -0.1325 -0.1325 -0.1325 0.0000 30.0000 -0.1325 -0.1325 -0.1325 -0.0000 2、迭代法求方程: 迭代法输出结果: >> f=inline('x^2-sin(x+0.15)'); >> [x0,k]=iterate(f,1.96) x0 = NaN k = 12 >> [x0,k]=iterate(f,0,2) x0 = -0.1494 k = 1 >> [x0,k]=iterate(f,0.2) x0 = 0.3234 k = 500 >> [x0,k]=iterate(f,0.3) x0 = 0.3234 k = 500 >> [x0,k]=iterate(f,0.001) x0 = 0.3234 k = 500 (3)、误差分析:由二分法和迭代法输出结果可知,利用二分法求出的结果中,可以清楚看出方程等于零时的解,其误差比迭代法小。 d、g(x)=x x-cos(x) 二分法求方程: (1)、在matlab的命令窗口中输入命令: >> fplot('[x^(x-cos(x))]',[-1,1]);grid 得下图: 由上图可得知:方程在[-1,1]区间有根。 (2)、二分法输出结果 >> f='x^(x-cos(x))' f = x^(x-cos(x)) >> bisect(f,-0.1,0.1,20,10^(-12)) 1.0000 -0.1000 0.1000 0 Inf 2.0000 -0.1000 0 -0.0500 -22.8740 + 3.5309i 3.0000 -0.0500 0 -0.0250 -43.6821 + 3.3947i 4.0000 -0.0250 0 -0.0125 -84.4110 + 3.2958i 1.0e+002 * 0.0500 -0.0001 0 -0.0001 -1.6511 + 0.0323i 1.0e+002 * 0.0600 -0.0001 0 -0.0000 -3.2580 + 0.0319i 1.0e+002 * 0.0700 -0.0000 0 -0.0000 -6.4648 + 0.0317i 1.0e+003 * 0.0080 -0.0000 0 -0.0000 -1.2872 + 0.0032i 1.0e+003 * 0.0090 -0.0000 0 -0.0000 -2.5679 + 0.0032i 1.0e+003 * 0.0100 -0.0000 0 -0.0000 -5.1285 + 0.0031i 1.0e+004 * 0.0011 -0.0000 0 -0.0000 -1.0249 + 0.0003i 1.0e+004 * 0.0012 -0.0000 0 -0.0000 -2.0490 + 0.0003i 1.0e+004 * 0.0013 -0.0000 0 -0.0000 -4.0971 + 0.0003i 1.0e+004 * 0.0014 -0.0000 0 -0.0000 -8.1931 + 0.0003i 1.0e+005 * 0.0001 -0.0000 0 -0.0000 -1.6385 + 0.0000i 1.0e+005 * 0.0002 -0.0000 0 -0.0000 -3.2769 + 0.0000i 1.0e+005 * 0.0002 -0.0000 0 -0.0000 -6.5537 + 0.0000i 1.0e+006 * 0.0000 -0.0000 0 -0.0000 -1.3107 + 0.0000i 1.0e+006 * 0.0000 -0.0000 0 -0.0000 -2.6215 + 0.0000i 1.0e+006 * 0.0000 -0.0000 0 -0.0000 -5.2429 + 0.0000i 2、迭代法求方程: 迭代法输出结果: >> f=inline('x^2-sin(x+0.15)'); x0 = 0.3234 k = 500 >> [x0,k]=iterate(f,0.01) x0 = 0.3234 k = 500 >> [x0,k]=iterate(f,0.81) x0 = 0.3234 k = 500 >> [x0,k]=iterate(f,0.61) x0 = 0.3234 k = 500 (3)、误差分析:由二分法和迭代法输出结果可知,利用二分法求出的结果中,可以清楚看出方程等于零时的解,其误差比迭代法小。 实验十七牛顿迭代法 【实验目的】 1.了解牛顿迭代法的基本概念。 2.了解牛顿迭代法的收敛性和收敛速度。 3.学习、掌握MATLAB软件的有关命令。 【实验内容】 用牛顿迭代法求方程的近似根,误差不超过。 3210 ++-=3 10- x x x 【实验准备】 1.牛顿迭代法原理 2.牛顿迭代法的几何解析 3.牛顿迭代法的收敛性 4.牛顿迭代法的收敛速度 5.迭代过程的加速 6.迭代的MATLAB命令 MATLAB中主要用for,while等控制流命令实现迭代。 【实验重点】 1.牛顿迭代法的算法实现 2.牛顿迭代法收敛性和收敛速度 【实验难点】 1.牛顿迭代法收敛性和收敛速度 【实验方法与步骤】 练习1用牛顿迭代法求方程在x=0.5附近的近似 3210 ++-= x x x 根,误差不超过。 310-牛顿迭代法的迭代函数为 322()1()()321 f x x x x g x x x f x x x ++-=-=-'++相应的MATLAB 代码为 >>clear; >>x=0.5; >>for i=1:3 >>x=x-(x^3+x^2+x-1)/(3*x^2+2*x+1) >>end 可算的迭代数列的前3项0.5455,0.5437,0.5437。经三次迭代,就大大超过了精度要求。 练习2 用牛顿迭代法求方程的近似正实根,由此建2(0)x a a =>立一种求平方根的计算方法。 由计算可知,迭代格式为,在实验12的练习4中1()()2a g x x x =+已经进行了讨论。 【练习与思考】 1.用牛顿迭代法求方程的近似根。 ln 1x x =2.为求出方程的根,在区间[1,2]内使用迭代函数进行310x x --=迭代,纪录迭代数据,问迭代是否收敛?对迭代进行加速,对比加速前的数据,比较加速效果。 3.使用在不动点的泰勒公式,证明牛顿迭代法收敛原理。*x 用二 4224min ()f t t t t =--[,.]t ∈内的极小值点,要求准 1. function [t d]=erfenfa(a,b) k=1; %记录循环次数 while abs(a-b)>0.0005 c=(a+b)/2; C(k)=c; %存储每次循环中点c 的值 if ff(c)<0 a=c; end if ff(c)==0 t1=c; break ; end if ff(c)>0 b=c; end k=k+1; end t=(a+b)/2; %最终符合要求的值 d=f(t); %最优解 C k function y=f(t) y=t^4-2*t^2-4*t; function y=ff(t) y=4*t^3-4*t-4; 运行结果 >> [t d]=erfenfa(1,1.5) C = Columns 1 through 9 1.2500 1.3750 1.3125 1.3438 1.3281 1.3203 1.3242 1.3262 1.3252 Column 10 1.3247 k = 11 t = 1.3250 d = -5.7290 2.黄金分割法 f (x)=x3-2x+1 初始区间[0, 3],收敛精度0.5 function [t,f]=huangjinfenge(a,b) m=1-(sqrt(5)-1)/2; t2=a+m*(b-a) f2=g(t2); t1=a+b-t2 f1=g(t1); while abs(t1-t2)>0.5 if f1 解线性方程组的迭代法 1.rs里查森迭代法求线性方程组Ax=b的解 function[x,n]=rs(A,b,x0,eps,M) if(nargin==3) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值elseif(nargin==4) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1; %迭代过程 while(tol>eps) x=(I-A)*x0+b; n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x; if(n>=M) disp('Warning:迭代次数太多,可能不收敛!'); return; end end 2.crs里查森参数迭代法求线性方程组Ax=b的解 function[x,n]=crs(A,b,x0,w,eps,M) if(nargin==4) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值 elseif(nargin==5) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1; %迭代过程 while(tol>eps) x=(I-w*A)*x0+w*b; n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x; if(n>=M) disp('Warning:迭代次数太多,可能不收敛!'); return; end end 3.grs里查森迭代法求线性方程组Ax=b的解 function[x,n]=grs(A,b,x0,W,eps,M) if(nargin==4) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值 elseif(nargin==5) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1;%前后两次迭代结果误差 %迭代过程 while(tol>eps) x=(I-W*A)*x0+W*b;%迭代公式 n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x; if(n>=M) disp('Warning:迭代次数太多,可能不收敛!'); return; end end 4.jacobi雅可比迭代法求线性方程组Ax=b的解 function[x,n]=jacobi(A,b,x0,eps,varargin) if nargin==3 eps=1.0e-6; M=200; elseif nargin<3 error return elseif nargin==5 M=varargin{1}; end D=diag(diag(A));%求A的对角矩阵 L=-tril(A,-1);%求A的下三角阵 function iteration A=[10,1,2,3,4; 1,9,-1,2,-3; 2,-1,7,3,-5; 3,2,3,12,-1; 4,-3,-5,-1,15]; b=[12,-27,14,-17,12]'; x0=[0,0,0,0,0]'; tol=1e-12; disp('jacobi迭代法的结果和次数如下:') [x,k]=Fjacobi(A,b,x0,tol) disp('G-S迭代法的结果和次数如下:':') [x,k]=Fgseid(A,b,x0,tol) disp('超松弛的结果和次数如下:':') [x,k]=Fsor(A,b,x0,1.2,tol) disp('共轭梯度法的结果和次数如下:':') [x,k]=Fcg(A,b,x0,tol) %jacobi迭代法 function [x,k]=Fjacobi(A,b,x0,tol) max=300; D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); B=D\(L+U); f=D\b; x=B*x0+f; k=1; while norm(x-x0)>=tol x0=x; x=B*x0+f; k=k+1; if(k>=max) disp('μü′ú3?1y300′?£?·?3ì×é?é?ü2?ê?á2'); return; end end %G-S迭代法 function [x,k]=Fgseid(A,b,x0,tol) max=300; D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); G=(D-L)\U; f=(D-L)\b; x=G*x0+f; k=1; while norm(x-x0)>=tol x0=x; x=G*x0+f; k=k+1; if(k>=max) disp('μü′ú3?1y300′?£?·?3ì×é?é?ü2?ê?á2'); return; end 要求: 下面分别使用雅克比迭代法和高斯-赛德尔迭代法求一个方程组的近似解用的线性方程组是按实验要求给的: 7*x1+x2+2*x3=10 x1+8*x2+2*x3=8 2*x1+2*x2+9*x3=6 雅克比迭代法的matlab代码:(老师写的) A=[7,1,2;1,8,2;2,2,9]; b=[10;8;6]; if(any(diag(A))==0) error('error,pause') end eps=input('误差限eps='); N=input('迭代次数N='); D=diag(diag(A)); B=inv(D)*(D-A); f=inv(D)*b; K=0; x0=zeros(size(b)); while 1 x1=B*x0+f K=K+1; fprintf('第-次迭代的近似解为',K) disp(x1'); if norm(x1-x0,inf) A=[7,1,2;1,8,2;2,2,9]; b=[10;8;6]; if(all(diag(A))==0) error('error,pause') end eps=input('误差限eps='); N=input('迭代次数N='); D=diag(diag(A)); B=inv(D)*(D-A); f=inv(D)*b; K=0; x0=zeros(size(b)); x00=x0; while 1 x11=B*x0+f; x00(1,1)=x11(1,1); x12=B*x00+f; x00(2,1)=x12(2,1); x13=B*x00+f; x00(3,1)=x13(3,1); x1=x00 K=K+1; fprintf('第-次迭代的近似解为',K) disp(x1'); if norm(x1-x0,inf) 实验一非线性方程的数值解法(一) 信息与计算科学金融崔振威201002034031一、实验目的: 熟悉二分法和简单迭代法的算法实现。 二、实验内容: 教材P40 2.1.5 三、实验要求 1根据实验内容编写二分法和简单迭代法的算法实现 2简单比较分析两种算法的误差 3试构造不同的迭代格式,分析比较其收敛性 (一)、二分法程序: function ef=bisect(fx,xa,xb ,n, delta) % fx是由方程转化的关于x的函数,有fx=0。 % xa解区间上限 % xb解区间下限 % n最多循环步数,防止死循环。 %delta为允许误差 x=xa;fa=eval(fx); x=xb;fb=eval(fx); disp(' [ n xa xb xc fc ]'); for i=1: n xc=(xa+xb)/2;x=xc;fc=eval(fx); X=[i,xa,xb,xc,fc]; disp(X), if fc*fa<0 xb=xc; else xa=xc; end if (xb-xa) k=0; while abs(x-xO)>eps & k 一、实验目的及题目 1.1 实验目的: (1)学会用高斯列主元消去法,LU 分解法,Jacobi 迭代法和Gauss-Seidel 迭代法解线性方程组。 (2)学会用Matlab 编写各种方法求解线性方程组的程序。 1.2 实验题目: 1. 用列主元消去法解方程组: 1241234 123412343421233234x x x x x x x x x x x x x x x ++=??+-+=??--+=-??-++-=? 2. 用LU 分解法解方程组,Ax b =其中 4824012242412120620266216A --?? ?- ?= ? ?-??,4422b ?? ? ?= ?- ?-?? 3. 分别用Jacobi 迭代法和Gauss-Seidel 迭代法求解方程组: 123234 1231234102118311210631125x x x x x x x x x x x x x -+=-??-+=-??-+=??-+-+ =? 二、实验原理、程序框图、程序代码等 2.1实验原理 2.1.1高斯列主元消去法的原理 Gauss 消去法的基本思想是一次用前面的方程消去后面的未知数,从而将方程组化为等价形式: 1111221122222n n n n nn n n b x b x b x g b x b x g b x g +++=??++=????= ? 这个过程就是消元,然后再回代就好了。具体过程如下: 对于1,2, ,1k n =-,若() 0,k kk a ≠依次计算 ()() (1)()()(1)()()/,,1, ,k k ik ik kk k k k ij ij ik kj k k k i i ik k m a a a a m a b b m b i j k n ++==-=-=+ 然后将其回代得到: ()() ()()()1/()/,1,2,,1 n n n n nn n k k k k k kj j kk j k x b a x b a x a k n n =+?=??=-=--? ? ∑ 以上是高斯消去。 但是高斯消去法在消元的过程中有可能会出现() 0k kk a =的情况,这时消元就无法进行了,即使主元数() 0,k kk a ≠但是很小时,其做除数,也会导致其他元素数量级的严重增长和舍入误差的扩散。因此,为了减少误差,每次消元选取系数矩阵的某列中绝对值最大的元素作为主元素。然后换行使之变到主元位置上,再进行销元计算。即高斯列主元消去法。 2.1.2直接三角分解法(LU 分解)的原理 先将矩阵A 直接分解为A LU =则求解方程组的问题就等价于求解两个三角形方程组。 直接利用矩阵乘法,得到矩阵的三角分解计算公式为: 1111111 11 1,1,2,,/,2,,,,,1,,,2,3, ()/,1,2, ,i i i i k kj kj km mj m k ik ik im mk kk m u a i n l a u i n u a l u j k k n k n l a l u u i k k n k n -=-===?? ==?? =-=+??=??=-=++≠?? ∑∑且 由上面的式子得到矩阵A 的LU 分解后,求解Ux=y 的计算公式为 11 111,2,3,/()/,1,2, ,1 i i i ij j j n n nn n i i ij j ii j i y b y b l y i n x y u x y u x u i n n -==+=??? =-=?? =??? =-=--?? ∑∑ 以上为LU 分解法。 二分法的matlab主程序 function [k,x,wuca,yx]=erfen(a,b,abtol) a(1)=a; b(1)=b; ya=fun(a(1)); yb=fun(b(1)); %程序中调用的fun.m 为函数if ya* yb>0, disp('注意:ya*yb>0,请重新调整区间端点a和b.'), return end max1=-1+ceil((log(b-a)- log(abtol))/ log(2)); % ceil是上取整for k=1: max1+1 a;ya=fun(a); b;yb=fun(b); x=(a+b)/2; yx=fun(x); wuca=abs(b-a)/2; k=k-1; [k,a,b,x,wuca,ya,yb,yx] if yx==0 a=x; b=x; elseif yb*yx>0 b=x;yb=yx; else a=x; ya=yx; end if b-a< abtol , return, end end k=max1; x; wuca; yx=fun(x); 区间二分法: 与对分查找法相同 1 区间二分法求出的仅仅是方程的一个单根,如果方程有重根或者多个根时,在做区间二分法时就会出现分叉,这样方程有几个根,就会产生几个实数序列,每一个实数序列的极限便是方程的一个根 2 通常用区间二分法为一些迭代法提供靠近x^*的初始选代值; 3 区间二分法的缺点是不能求方程的复数根。 format long a=5; b=6; x1=a; x2=b; f1=4*cos(x1)+4*sin(x1)+0.5*x1-2; f2=4*cos(x2)+4*sin(x2)+0.5*x2-2; step=0.000001; ii=0; while abs(x1-x2)>step ii=ii+1; x3=(x1+x2)/2; f3=4*cos(x3)+4*sin(x3)+0.5*x3-2; if f3~=0 if f1*f3<0 x2=x3; else x1=x3; end end end x3 f=[4*cos(x3)+4*sin(x3)+0.5*x3] disp(['迭代次数:',num2str(ii),'次']) 牛顿迭代法求解: 在方程f(x)=0有实数根的情况下,若能够将方程等价地转化成x=g(x)的形式,然后取一个初始值x0代入x=g(x)的右端,算得x1=g(x0),再计算x2=g(x1),这样依次类推 x(k+1)=g(x(k)) 可以得到一个序列xk,通常称g(x)为迭代函数,序列xk为由迭代函数产生得迭代序列,x0为迭代初始值。 同一个方程,不同等价形式的转换产生的迭代法可能收敛,也有可能发散.关于迭代法的敛散性判定有下面的定理(也称李普希兹(Lipschitz定理): 如果迭代函数g(x)在区间[a,b]上连续,且满足以下条件, 1 对于任意的x=[a,b],有g(x)=[a,b] 1.编写程序:计算1/3+2/5+3/7+……+10/21 法一: s=0; for i=1:10 s=s+i/(2*i+1); end s s = 4.4096 法二: sum((1:10)./(3:2:21)) ans = 4.4096 2.编写程序:计算1~100中即能被3整除,又能被7整除的所有数之和。 s=0; for i=1:100 if mod(i,3)==0&&mod(i,7)==0 s=s+i; end,end s s = 210 3.画出y=n!的图(1<=n<=10),阶乘的函数自己编写,禁用MATLAB自带的阶乘函数。 x=1:10; for i=1:10 try y(i)=y(i-1)*i; catch y(i)=1; end,end plot(x,y) 12345678910 0.511.522.533.54x 10 6 4.一个数恰好等于它的因子之和,这个数就称为完数。例如,6的因子为1,2,3,而6=1+2+3,因此6就是一个完数。编程找出2000以内的所有完数。 g=[]; for n=2:2000 s=0; for r=1:n-1 if mod(n,r)==0 s=s+r; end end if s==n g=[g n]; end end g g =6 28 496 5.编写一个函数,模拟numel函数的功能,函数中调用size函数。 function y=numelnumel(x) m=size(x); y=m(1)*m(2); numelnumel([1 2 3;4 5 6]) ans = 6 6. 编写一个函数,模拟length函数的功能,函数中调用size函数。 function y=lengthlength(x) m=size(x); y=max(m(1),m(2)); lengthlength([1 2 3;4 5 6]) ans = 3 7.求矩阵rand(5)的所有元素和及各行平均值,各列平均值。 s=rand(5); sum=sum(sum(s)) mean2=mean(s,2) mean1=mean(s) sum = 13.8469 迭代法最佳松弛因子的选取 一、问题提出: 针对矩阵430341014A ?? ??=-?? ??-?? ,b=[24;30;-24],用SOR 迭代求解。并选出最佳松弛 因子。理论分析 1.24ω==≈。做出()L ωρ关于ω函数 的图像。 二、理论基础 选取分裂矩阵M 为带参数的下三角矩阵)(1 wL D w M -=, 其中w>0为可选择的松弛因子. 于是,由 ?????+=+f Bx x x k k ) ()1()0() (初始向量 (k=0,1,…,)可构造一个迭代法,其迭代矩阵为A wL D w I L w 1)(---≡ =).)1(()(1wU D w wL D +--- 从而得到解Ax=b 的主次逐次超松弛迭代法. 解Ax=b 的SOR 方法为 ?????+=+f Bx x x k k ) ()1()0() (初始向量 (k=0,1,…,) (1) 其中 w L =).)1(()(1wU D w wL D +---(2) b wL D w f 1)(--= 下面给出解Ax=b 的SOR 迭代法的分量计算公式.记 ,),...,,...,() () () (1)(T k n k i k k x x x x = 由(1)式可得 ,))1(()()()1(wb x wU D w x wL D k k ++-==-+ ).()()()1()()1(k k k k k Dx Ux Lx b w Dx Dx -+++=++ (3) 由此,得到解Ax=b 的SOR 方法的计算公式 ?????????==--+==∑∑-==++.),1,0;,...,2,1(/)(,),...,(11) (1)()1()0()0(1)0(为松弛因子 w k n i a x a x a b w x x x x x ii i j n i j k j ij k j ij i k i k i T n (4) 或 ?? ?? ? ??????==--=??+==∑∑-==++.,...),1,0;,...,2,1()/(,,),...,(.11)()1() () 1()0()0(1)0(为松弛因子w k n i a x a x a b w x x x x x x x i j n i j ii k j ij k j ij i i i k i k i T n (5) ※ 若要求选取出最佳松弛因子,则有两种方法: ⑴、 给出w 的最佳范围,当取不同的w 值时,会求出不同的谱半径R 的值, 然后判断出值最小的谱半径。那么这个最小的谱半径所对应的w ,即为所求最佳松弛因子。 ⑵、 给出w 的最佳范围,当取不同的w 值时,由(2)式进行迭代,看它们在 相同精度范围内的迭代次数,找出迭代次数最低的那一个,其所应用的w 即为最佳松弛因子。 三、实验内容: 从表格中可以看出,迭代次数随着松弛因子的增长而呈现先减后增的趋势,当谱半径最小时,其迭代次数最小。则表示出谱半径最小时,其松弛因子为最佳松弛因子。 matlab 迭代法[精品] 1. 矩阵 122,211,,,,,,,,,A,111A,222, 11,,,,,,,,221,,112,,,, 证明:求解以为系数矩阵线性方程组的Jacobi迭代式收敛的,而A1 Gauss-Seidel方法是发散的;求解以为系数矩阵线性方程组的A2实验名称Gauss-Seidel是收敛的,而Jacobi方法是发散的. 2. 矩阵 1aa,,,,Aaa,1 ,,,,aa1,, (a) 参数取什么值时,矩阵是正定的. a (b) 取什么值时,求以为系数矩阵线性方程组的Jacobi迭代式收aa 敛的. 1、根据迭代收敛性的充分必要条件来判断Jacobi迭代式与Gauss-Seide 迭代式的收敛性,迭代收敛性仅与方程组系数矩阵有关,与右端无关;而且不依赖于初值的选取。实验目的 2、根据矩阵的判断定理求得矩阵元素a的取值,同时根据矩阵线性方程组的Jacobi迭代式收敛的充分条件(严格对角占优)来求a得取值。 1、(1)检验线性方程组的Jacobi迭代式的收敛性: function jacobi(A) D=zeros(3); for i=1:3 D(i,i)=A(i,i); 实验内容end (算法、程B=D^(-1)*(D-A); 序、步骤和k=max(abs(eig(B))) 方法) if k<1 '该线性方程组的Jacobi迭代式是收敛的' else k>=1 '该线性方程组的Jacobi迭代式是发散的' end (2)检验线性方程组的Gauss-Seide迭代式的收敛性: function Gauss(A) D=zeros(3); L=zeros(3); U=zeros(3); for i=1:3 D(i,i)=A(i,i); end L(2:3,1)=A(2:3,1); L(3,2)=A(3,2); U(1,2:3)=A(1,2:3); U(2,3)=A(2,3); B=-(D+L)^(-1)*U; k=max(abs(eig(B))) if k<1 '该线性方程组的Gauss-Seidel迭代式是收敛的' else k>=1 '该线性方程组的Gauss-Seidel迭代式是发散的' end 2、(1)参数取什么值时,矩阵是正定的.(矩阵的特征值全为正) a >> syms a >> A=[1 a a;a 1 a;a a 1]; >> eig(A) ans = 2*a+1 1-a 实验四 姓名:木拉丁。尼则木丁班级:信计08-2 学号:20080803405 实验地点:新大机房 实验目的:通过本实验学习利用MATLAB不动点迭代法,抛物线法,斯特芬森迭代法解非线性方程组,及其编程实现,培养编程与上机调试能力。 实验要求:①上机前充分准备,复习有关内容,写出计算步骤,查对程序; ②完成实验后写出完整的实验报告,内容应该包括:所用的算法语言, 算法步骤陈述,变量说明,程序清单,输出计算结果,结果分析等等; ③用编好的程序在Matlab环境中执行。 迭代法 MATLAB程序: function pwxff(f,x0,x1,x2,d,n) f=inline(f); x(1)=x0; x(2)=x1; x(3)=x2; w1=(f(x(2))-f(x(3)))/(x(2)-x(3)); t1=(f(x(1))-f(x(3)))/(x(1)-x(3)); t2=(f(x(1))-f(x(2)))/(x(1)-x(2)); w2=1/(x(1)-x(2))*(t1-t2); w=w1+w2*(x(3)-x(2)); for k=3:n x(k+1)=x(k)-2*f(x(k))/(w+sqrt(w^2-4*f(x(k))*w2)); if abs(x(k+1)-x(k)) 实验一非线性方程的数值解法(一)信息与计算科学金融崔振威201002034031 一、实验目的: 熟悉二分法和简单迭代法的算法实现。 二、实验内容: 教材P40 2.1.5 三、实验要求 1 根据实验内容编写二分法和简单迭代法的算法实现 2 简单比较分析两种算法的误差 3 试构造不同的迭代格式,分析比较其收敛性 (一)、二分法程序: function ef=bisect(fx,xa,xb,n,delta) % fx是由方程转化的关于x的函数,有fx=0。 % xa 解区间上限 % xb 解区间下限 % n 最多循环步数,防止死循环。 %delta 为允许误差 x=xa;fa=eval(fx); x=xb;fb=eval(fx); disp(' [ n xa xb xc fc ]'); for i=1:n xc=(xa+xb)/2;x=xc;fc=eval(fx); X=[i,xa,xb,xc,fc]; disp(X), if fc*fa<0 xb=xc; else xa=xc; end if (xb-xa) k=0; while abs(x-x0)>eps & k 第六章线性方程组的迭代解法 2015年12月27日17:12 迭代法是目前求解大规模稀疏线性方程组的主要方法之一。包括定常迭代法和不定常迭代法,定常迭代法的迭代矩阵通常保持不变,包括有雅可比迭代法(Jacobi)、高斯-塞德尔迭代法(Gauss-Seidel)、超松弛迭代法(SOR) 1.雅可比迭代法(Jacobi) A表示线性方程组的系数矩阵,D表示A的主对角部分,L表示下三角部分,U表示上三角部分。 A=D+L+U 要解的方程变为Dx+Lx+Ux=b x=D^(-1)(b-(L+U)x) 所以Jocabi方法如下: Matlab程序 function [x,iter] =jacobi(A,b,tol) D=diag(diag(A)); L=D-tril(A); U=D-triu(A); x=zeros(size(b)); for iter=1:500 x=D\(b+L*x+U*x); error=norm(b-A*x)/norm(b); if(error matlab 迭代法代码 1、%用不动点迭代法求方程x-e A x+4=0的正根与负根,误差限是 10A-6% disp(' 不动点迭代法 '); n0=100; p0=-5; for i=1:n0 p=exp(p0)-4; if abs(p-p0)<=10(6) if p<0 disp('|p-p0|=') disp(abs(p-p0)) disp(' 不动点迭代法求得方程的负根为 :') disp(p); break; else disp(' 不动点迭代法无法求出方程的负根 .') end else p0=p; end end if i==n0 disp(n0) disp(' 次不动点迭代后无法求出方程的负根') end p1=1.7; for i=1:n0 pp=exp(p1)-4; if abs(pp-p1)<=10(6) if pp>0 disp('|p-p1|=') disp(abs(pp-p1)) disp(' 用不动点迭代法求得方程的正根为 ') disp(pp); else disp(' 用不动点迭代法无法求出方程的正根 '); end break; else p1=pp; end end if i==n0 disp(n0) disp(' 次不动点迭代后无法求出方程的正根 ') end 2、%用牛顿法求方程x-e A x+4=0的正根与负根,误差限是disp(' 牛顿法') n0=80; p0=1; for i=1:n0 p=p0-(p0-exp(p0)+4)/(1-exp(p0)); if abs(p-p0)<=10(6) disp('|p-p0|=') disp(abs(p-p0)) disp(' 用牛顿法求得方程的正根为 ') disp(p); break; else p0=p; end end if i==n0 disp(n0) disp(' 次牛顿迭代后无法求出方程的解 p1=-3; for i=1:n0 p=p1-(p1-exp(p1)+4)/(1-exp(p1)); 10A-6 ') end 1.列主元高斯消去法 M文件 function[x]=gauss(a,b) n=length(a); x=zeros(n,1); a=[a b]; for k=1:n-1 max=k; for i=k+1:n if a(i,k)>a(max,k) max=i; end end temp=a(k,k:n+1); a(k,k:n+1)=a(max,k:n+1); a(max,k:n+1)=temp; for i=k+1:n a(i,k)=-a(i,k)/a(k,k); a(i,k+1:n+1)=a(i,k+1:n+1)+a(i,k)*a(k,k+1:n+1); end end x(n,1)=a(n,n+1)/a(n,n); for i=n-1:-1:1 sum=0; for j=i+1:n sum=sum+x(j,1)*a(i,j); end x(i,1)=(a(i,n+1)-sum)/a(i,i); end Matlab运行结果 2.LU三角分解法 M文件 function y=LU(A,B); n=length(A); A=[A B]; for k=1:n-1; for i=k:n; if(abs(A(i,k))==max(abs(A(k:n,k)))) P(k)=i; temp=A(k,:); A(k,:)=A(i,:); A(i,:)=temp; end end for j=k+1:n; A(j,k)=A(j,k)/A(k,k); A(j,k+1:n+1)=A(j,k+1:n+1)-A(j,k)*A(k,k+1:n+1); end end P(n)=n; L(1,1)=1; L(2:n,1)=A(2:n,1); L(1,2:n)=0; U(1,1)=A(1,1); U(2:n,1)=0; U(1,2:n)=A(1,2:n); for i=2:n; L(i,1:i-1)=A(i,1:i-1); L(i,i)=1; L(i,i+1:n)=0; U(i,1:i-1)=0; U(i,i:n)=A(i,i:n); end x(n) = A(n,n+1)/U(n,n); for k = n-1:-1:1 x(k)=A(k,n+1); for p=n:-1:k+1; 二分法、简单迭代法的m a t l a b代码实现 实验一非线性方程的数值解法(一)信息与计算科学金融崔振威 201002034031 一、实验目的: 熟悉二分法和简单迭代法的算法实现。 二、实验内容: 教材P40 2.1.5 三、实验要求 1 根据实验内容编写二分法和简单迭代法的算法实现 2 简单比较分析两种算法的误差 3 试构造不同的迭代格式,分析比较其收敛性 (一)、二分法程序: function ef=bisect(fx,xa,xb,n,delta) % fx是由方程转化的关于x的函数,有fx=0。 % xa 解区间上限 % xb 解区间下限 % n 最多循环步数,防止死循环。 %delta 为允许误差 x=xa;fa=eval(fx); x=xb;fb=eval(fx); disp(' [ n xa xb xc fc ]'); for i=1:n xc=(xa+xb)/2;x=xc;fc=eval(fx); X=[i,xa,xb,xc,fc]; disp(X), if fc*fa<0 xb=xc; else xa=xc; end if (xb-xa) 二分法 二分法基本思路 一般地,对于函数f(x),如果存在实数C,当x=c时,若f(c)=O,那么把x=c叫做函数f(x)的零点。解方程即要求f(x)的所有零点。 假定f(x)在区间(x,y)上连续 先找到a、b属于区间(x, y),使f(a),f(b)异号,说明在区间(a,b)内一定有零点,然后求 f[(a+b)/2], 现在假设f(a) 解线性方程组b AX =的迭代法是从初始解出发,根据设计好的步骤用逐次求出的近似解逼近精确解.在第三章中介绍的解线性方程组的直接方法一般适合于A 为低阶稠密矩阵(指n 不大且元多为非零)的情况,而在工程技术和科学计算中常会遇到大型稀疏矩阵(指n 很大且零元较多)的方程组,迭代法在计算和存贮两方面都适合后一种情况.由于迭代法是通过逐次迭代来逼近方程组的解,所以收敛性和收敛速度是构造迭代法时应该注意的问题.另外,因为不同的系数矩阵具有不同的性态,所以大多数迭代方法都具有一定的适用范围.有时,某种方法对于一类方程组迭代收敛,而对另一类方程组迭代时就发散.因此,我们应该学会针对具有不同性质的线性方程组构造不同的迭代. 4.1 迭代法和敛散性及其MATLAB 程序 4.1.2 迭代法敛散性的判别及其MATLAB 程序 根据定理4.1和谱半径定义,现提供一个名为pddpb.m 的M 文件,用于判别迭代公H=eig(B);mH=norm(H,inf); if mH>=1 disp('请注意:因为谱半径不小于1,所以迭代序列发散,谱半径mH 和B 的所 有的特征值H 如下:') else disp('请注意:因为谱半径小于1,所以迭代序列收敛,谱半径mH 和B 的所有 的特征值H 如下:') end mH 4.1.3 与迭代法有关的MATLAB 命令 (一) 提取(产生)对角矩阵和特征值 可以用表4–1的MATLAB 命令提取对角矩阵和特征值. (二) 提取(产生)上(下)三角形矩阵 可以用表4–2的MATLAB命令提取矩阵的上三角形矩阵和下三角形矩阵. (三)稀疏矩阵的处理 对稀疏矩阵在存贮和运算上的特殊处理,是MA TLAB进行大规模科学计算时的特点和优势之一.可以用表4–3的MATLAB命令,输入稀疏矩阵的非零元(零元不必输入),即可进行运算. 4.2 雅可比(Jacobi)迭代及其MATLAB程序 4.2.2 雅可比迭代的收敛性及其MATLAB程序 [n m]=size(A); for j=1:m a(j)=sum(abs(A(:,j)))-2*(abs(A(j,j))); end for i=1:n if a(i)>=0 disp('请注意:系数矩阵A不是严格对角占优的,此雅可比迭代不一定收敛') return end end if a(i)<0 disp('请注意:系数矩阵A是严格对角占优的,此方程组有唯一解,且雅可比迭代收敛') end 例4.2.2 用判别雅可比迭代收敛性的MATLAB主程序,判别由下列方程组的雅可比迭matlab实验十七__牛顿迭代法(可打印修改)
matlab编程实现二分法,牛顿法,黄金分割法,最速下降matlab程序代码
MATLAB代码 解线性方程组的迭代法
jacobi G-S,超松弛迭代法MATLAB程序
MATLAB样例之雅克比迭代法
二分法、简单迭代法的matlab代码实现
lu分解法、列主元高斯法、jacobi迭代法、gaussseidel法的原理及matlab程序
matlab 二分法汇总
matlab程序设计例题及答案
MATLAB实现迭代法最佳松弛因子的选取
matlab 迭代法[精品]
【良心出品】不动点迭代法matlab程序
二分法、简单迭代法的matlab代码实现教学文案
线性方程组的迭代解法(Matlab)
matlab迭代法代码
方程组的各种解法法的Matlab程序及运行结果
二分法、简单迭代法的matlab代码实现教学文案
二分法matlab程序
迭代解法的matlab实现