matlab(计算方法的MATLAB实现)

5.4000 7.0500 7.5000
• >> x=[1 2 4 6 8 9 10 13 15 16]; • >> y=[5 7 8 10 13 14 15 17 19 20]; • >> x1=[1.2 2.1 3]; • >> y1=interp1(x,y,x1,'linear')
• y1 = •
• Gauss-Saidel迭代法 • 编写gauss函数，并使之计算 • >> a=[10 -2 -1;-2 10 -1;-1 -2 5]; • >> b=[3 15 10]'; • >> s=gauss(a,b,[0 0 0]',eps) •s= • 1 • 2 • 3
• • • • • • • • • • • • • • • • • • • • •
• 2、线性方程组求解中的变换
• 上三角变换 • U=triu(x)返回矩阵x的上三角部分； • U=triu(x,k)返回第k条对角线以上部分的元
素。
• >> a=[12 -3 3;-18 3 -1;1 1 1]; • >> triu(a) • ans =
• • •
12 0 0 -3 3 3 -1 0 1
MATLAB 7.0从入 门到精通
主要讲述内容
• 第1章 MATLAB简介 • 第2章 数值运算 • 第3章 单元数组和结构 • 第4章 字符串 • 第5章 符号运算 • 第6章 MATLAB绘图基础 • 第7章 程序设计 • 第8章 计算方法的MATLAB实现 • 第9章 优化设计 • 第10章 Simulink仿真初探
function s=gauss(a,b,x0,eps) %gauss-seidel迭代法皆线性方程组 %a为系数矩阵，b为方程组ax=b中的右边的矩阵b，x0为迭代初值 if nargin==3 eps=1.0e-6; elseif nargin<3 error return end L=tril(a,-1);%求出严格下三角矩阵 D=diag(diag(a));%求出对角矩阵 U=triu(a,1);%求出严格上三角矩阵 C=inv(D+L); B=-C*U; f=C*b; s=B*x0+f; while norm(s-x0)>=eps x0=s; s=B*x0+f; end return
• 8.1 方程求根
• >> fzero('x^3-3*x-1',2) • ans = • 1.8794 • >> fzero('x^3-3*x-1',[2,4]) • ??? Error using ==> fzero • The function values at the interval
endpoints must differ in sign. • >> fzero('x^3-3*x-1',[1,4]) • ans = • 1.8794
• • •
12 -18 1
-3 3 1
0 -1 1
• >> tril(a,-1) • ans =
• • •
0 -18 1
0 0 1
0 0 0
• 3、迭代解法
• 迭代解法非常适用于求解大型稀疏系数矩阵的
方程组，在线性方程组常用的迭代解法主要有 Jacobi迭代法、Gauss-Seidel迭代法。
• Jacobi迭代法 • 编写jacobi函数，并使之计算 • >> a=[10 -2 -1;-2 10 -1;-1 -2 5]; • >> b=[3 15 10]'; • >> s=jacobi(a,b,[0 0 0]',eps) •s= • 1 • 2 • 3
• Method可取如下的值： • ‘nearest’最近插值 • ‘linear’线性插值 • ‘spline’三次样条插值 • ‘cubic’三次插值 • Method默认值为线性插值，上述插值要求
向量x单调。
• >> x=[1 2 4 6 8 9 10 13 15 16]; • >> y=[5 7 8 10 13 14 15 17 19 20]; • >> x1=[1.2 2.1 3]; • >> y1=interp1(x,y,x1) • y1 = •
• ans = •
1 1
• • • • • • • • • • • • • • • • • • • • •
function s=newtoniterate(x,eps) %newton迭代法解非线性方程组，x为迭代初值，eps为允许误差 if nargin==1 eps=1.0e-6; elseif nargin<1 error return end x1=fx1(x);%非线性方程组 x2=-dfx1(x);%非线性方程组导数 x3=inv(x2); x0=x3*x1'; while norm(x0)>=eps x=x0'+x; x1=fx1(x); x2=-dfx1(x); x3=inv(x2); x0=x3*x1'; end s=x0'+x; return
• >> triu(a,1)
• ans =
• • •
0 0 0 -3 0 0 3 -1 0
• >> triu(a,-1)
• ans = • • •
12 -18 0 -3 3 1 3 -1 1
• 对角变换 • U=diag(x)返回矩阵x主对角线上的元素，返
回结果是一列向量形式； • U=diag(x,k)返回第k条对角线上的元素值。 • 当x为向量时生成矩阵。
• Newton迭代法 • 用newton迭代法求解下面的方程组
x1 10 x1 x2 8 0 2 x1 x2 x1 10 x2 8 0
2 2
• 首先，编写上述非线性方程组的M文件
fx1.m • 然后，编写上述非线性方程组导数的M文件 dfx1.m • >> newtoniterate([0 0])
must differ in sign.
• 1、直接解法
• 8.2 线性方程组数值解法
• A*X=B • X=A\B或x=inv(A)*B • 例： • 12x(1)-3x(2)+3x(3)=15 • -18x(1)+3x(2)-x(3)=-15 • X(1)+x(2)+x(3)=6
• >> a=[12 -3 3;-18 3 -1;1 1 1]; • >> b=[15;-15;6]; • >> x=a\b •x= • 1.0000 • 2.0000 • 3.0000 • >> inv(a)*b • ans = • 1 • 2 • 3
x1 10 x1 x2 8 0 2 x1 x2 x1 10 x2 8 0
2 2
• 首先编写上述非线性方程组的M文件fx.m • >> staticiterate([0 0])
• ans =
•
1.0000 1.0000
• function s=staticiterate(x,eps) • %不动点迭代法解非线性方程组，x为迭代初值，eps为允 • • • • • • • • • • • • •
• >> a=[12 -3 3;-18 3 -1;1 1 1]; • >> diag(a)
• ans =
• • •
12 3 1
• >> a=[12 -3 3;-18 3 -1;1 1 1]; • >> diag(a,1)
• ans =
• •
-3 -1
• >> a=[12 -3 3;-18 3 -1;1 1 1]; • >> diag(a,-1) • ans = • •
-18 1
• 下三角变换 • U=tril(x)返回矩阵x的下三角部分； • U=tril(x,k)返回第k条对角线以上下部分的元
素。
• >> a=[12 -3 3;-18 3 -1;1 1 1]; • >> tril(a)
• ans =
• • •
12 -18 1
Biblioteka Baidu
0 3 1
0 0 1
• >> tril(a,1) • ans =
第8章 计算方法的MATLAB实现
• roots见多项式求根；roots(多项式系数矩阵) • fzero可求解非线性方程 • fzero(‘function’,x0)其中function为求解的方程，
x0为估计的根，x0可为标量或长度为2的向量， 为向量时函数的两端的值应该符号相反，此 时求区间上的解。只能求解x0附近的一个解。 即使在某个区间内有多个解，但是区间端点 符号相同的话仍然出错。
许误差 if nargin==1 eps=1.0e-6; elseif nargin<1 error return end xx=fx(x);%第一次迭代 while norm(xx-x)>=eps x=xx; xx=fx(x); end s=xx; return
• function y=fx(x) • y(1)=0.1*(x(1)*x(1)+x(2)*x(2)+8); • %y(1)=x(1)*x(1)+x(2)*x(2)-10*x(1)+8; • y(2)=0.1*(x(1)*x(2)*x(2)+x(1)+8); • %y(2)=x(1)*x(2)*x(2)+x(1)-10*x(2)+8; • y=[y(1) y(2)];
• 1、一维线性插值
• 8.4 插值与拟合
• yi=interp1(x,y,xi)返回在插值向量xi处的函数
向量yi，它是根据向量x和y插值而来。如y是矩 阵，则对y每一列进行插值，如xi中元素不在x 内，返回NaN。 • yi=interp1(y,xi)省略x，表示x=1:N，此时N为 向量y的长度或为矩阵y的行数。 • yi=interp1(x,y,xi,’method’)表示用method指定 的插值方法进行插值。
5.4000 7.0500 7.5000
• >> x=[1 2 4 6 8 9 10 13 15 16]; • >> y=[5 7 8 10 13 14 15 17 19 20]; • >> x1=[1.2 2.1 3]; • >> y1=interp1(x,y,x1,'nearest')
• 8.3 非线性方程组数值解法
• 与线性方程组的求解一样，非线性方程组的
求解也是应用很广泛的课题。一般情况，非 线性方程组的数据值解法主要采用迭代法来 求解。比较常用的迭代法主要有不动点迭代 法、Newton迭代法、拟Newton迭代法等几 种方法。
• 不动点迭代 • 编写不动点迭代函数，并使之计算 • 用不动点迭代法求解下面的方程组
• function y=fx1(x) • y(1)=x(1)*x(1)-10*x(1)+x(2)*x(2)+8; • y(2)=x(1)*x(2)*x(2)-10*x(2)+x(1)+8; • y=[y(1) y(2)];
• function y=dfx1(x) • y(1)=2*x(1)-10; • y(2)=2*x(2); • y(3)=x(2)*x(2)+1; • y(4)=2*x(1)*x(2)-10; • y=[y(1) y(2);y(3) y(4)];
• • • • • • • • • • • • • • • • • • • • •
function s=jacobi(a,b,x0,eps) %jacobi迭代法皆线性方程组 %a为系数矩阵，b为方程组ax=b中的右边的矩阵b，x0为迭代初值 if nargin==3 eps=1.0e-6; elseif nargin<3 error return end D=diag(diag(a));%求出对角矩阵 D=inv(D);%求出对角矩阵的逆 L=tril(a,-1);%求出严格下三角矩阵 U=triu(a,1);%求出严格上三角矩阵 B=-D*(L+U); f=D*b; s=B*x0+f; while norm(s-x0)>=eps x0=s; s=B*x0+f; end return
• >> fzero('x^2-3*x+2',0) • ans = • 1.0000 • >> fzero('x^2-3*x+2',3) • ans = • 2.0000 • >> fzero('x^2-3*x+2',[0,4]) • ??? Error using ==> fzero • The function values at the interval endpoints