现代数值计算方法(MATLAB版)第3章(2)
第3章MATLAB计算
1、一维插值运算
例子：某观测站测得某日6:00时至18:00时之间每隔两小时的室 外温度（℃）如表所示，用3次样条插值法分别求得该日室外 6:30时至17:30时之间每隔两小时各点的近似温度（℃）。
时间h 室内温度t1 6 18.0 8 20.0 10 22.0 12 25.0 14 30.0 16 28.0 18 24.0
Y=T*A
1、 回归法曲线拟合
y a0 a1et a2tet
>> T=[ones(size(t)), exp(-t),t.*exp(-t)] T= 1.0000 1.0000 0 1.0000 0.7408 0.2222 1.0000 0.4493 0.3595 1.0000 0.3329 0.3662 1.0000 0.2019 0.3230 1.0000 0.1003 0.2306 >> A=X\Y A= 1.3974 -0.8988 0.4097
例6-27
>> t=[0 0.3 0.8 1.1 1.6 2.3]'; >> y=[0.5 0.82 1.14 1.25 1.35 1.4]';
>> plot(t,y,'r*')
>> grid on
1、 回归法曲线拟合
y a0 a1t a2t 2
>> T=[ones(size(t)), t,t.^2] T = 1.0000 0 0 1.0000 0.3000 0.0900 1.0000 0.8000 0.6400 1.0000 1.1000 1.2100 1.0000 1.6000 2.5600 1.0000 2.3000 5.2900 >> A=T\Y A = 0.5318 0.9191 -0.2387
matlab学习3-数值计算
六、矩阵元素之间的逻辑运算
一、矩阵的构造
1、向量的构造
向量是1×N（ N×1 ）的特殊矩阵，称为N维向量。
是一种特殊的矩阵 （1）逐个输入法：x=[ ] 行向量：数据元素之间均用空格（或逗号）隔开； 例:x1=[2 3 sqrt(3) 5] 列向量：数据元素之间均用分号隔开 例:x2=[2;3;sqrt(3);5] 注：行向量和列向量之间的转换“ ’ ”
第二章
基本数值计算
第一节 简单的数学运算
第二节 MATLAB数值计算基础
第三节 MATLAB数值分析与多项式计算
第一节 简单的数学运算 一、常用的数学运算符 二、Matlab 语言规则 三、常用操作命令和键盘技巧 四、常量和变量 五、函数
一、常用的数学运算符
1、Matlab 的数学运算定义在复数域上。
example3
2、矩阵的基本运算： （1）标量与矩阵的数运算和数学函数对矩阵的运算等 于对矩阵的每一个元素的运算。 a=[1 2 3];b=a+100 b= 101 102 103 （2）进行矩阵加减时，参与运算的矩阵必须同维。 （3）进行矩阵乘法时， A的行数=B列数。 左乘与右乘不同：一般A*B不等于B*A 若A*B等于B*A，则称A，B对易 （4）幂运算A^n
2、对矩阵（A）的部分操作：
函数
Fliplr(A)
功能
矩阵左右翻转
函数
Tiag(A,k)
功能
取矩阵对角线 元素
Flipud(A)
Flipdim(A, m) Rot(A,k)
矩阵上下翻转
矩阵沿特定 维（m）翻转 矩阵逆时针旋 转k*90度
Tril(A,k)
Triu(A,k)
取矩阵的下三 角部分
第3章 MATLAB数值计算
• 2．利用矩阵三角分解LU、正交分解QR和特 征值分解求线性方程组的解 • 矩阵的三角分解LU分解、正交分解QR分解 和特征值分解已经在第2章讲过，在此不再 细讲。例如经过LU分解，方程组AX=B的系 数矩阵A可以化为A=LU，则方程组的解为： X= U-1（L-1B）
2013年7月27日星期六
>> clear >> a=magic(6) a=35 1 6 26 19 24 3 32 7 21 23 25 31 9 2 22 27 20 8 28 33 17 10 15 30 5 34 12 14 16 4 36 29 13 18 11 >> max(a,[],2) %求每行最大元素 ans=35 32 31 33 34 36
3.4 数据的数值计算
• MATLAB可以通过调用函数快速实现数据的 统计与分析、求解向量的内积与正交、数 据的插值与拟合以及插值拟合曲线的绘制 。
2013年7月27日星期六
31
3.4.1 数据统计与分析
• •
•
1．求矩阵最大元素max()和最小元素min() 【例3.13】求矩阵a=magic(6)的每行及每列的最大和最小元素，并求整个矩阵 的最大元素和最小元素。 解 在MATLAB命令提示符下输入：
2013年7月27日星期六
29
• 例：求有理分式f(x)=2x3+5x+3,g(x)=x+4的导 函数 • >>f=[2 0 5 3]; • df=polyder(f) • g=[1,4]; • dp=polyder(f,g) • [p,q]=polyder(f,g)
2013年7月27日星期六
30
3.3.2 多项式求根
matlab第3章 matlab3-数值计算
4. 矩阵的其它运算
inv —— 矩阵求逆 det —— 行列式的值 eig —— 矩阵的特征值 diag —— 对角矩阵
’ —— 矩阵转置 sqrt —— 矩阵开方
关系运算
关系符号 < <= > >= == ~= 意义 小于 小于或等于 大于 大于或等于 等于 不等于
5. 矩阵的数组运算(点运算)
d=[-1;0;2];f=pi*d f = -3.1416 0 6.2832 矩阵除的运算在线性代数中没有， 有矩阵逆的运算，在matlab中有两 种矩阵除运算
3. 矩阵乘方—— a^n,a^p,p^a
a ^ p —— a 自乘p次幂
方阵 >1的整数 的整数
如果a是矩阵，p是标量,a^p使用特征值 和特征向量自乘到p次幂；如a,p都是矩 阵，a^p则无意义。
3. 矩阵的修改
直接修改 可用↑键找到所要修改的矩阵，用←键 移动到要修改的矩阵元素上即可修改。 指令修改 可以用A(∗,∗)= ∗ 来修改。
例如 a=[1 2 0;3 0 5;7 8 9] a =1 2 0 3 0 5 7 8 9 a(3,3)=0 a =1 2 0 3 0 5 7 8 0
二、矩阵运算
1. 矩阵加、减（＋,－）运算
规则： 相加、减的两矩阵必须有相同的行和 列两矩阵对应元素相加减。 允许参与运算的两矩阵之一是标量。 标量与矩阵的所有元素分别进行加 减操作。
2. 矩阵乘（∗）运算
规则： A矩阵的列数必须等于B矩阵的行数 标量可与任何矩阵相乘。 a=[1 2 3;4 5 6;7 8 0];b=[1;2;3];c=a*b c =14 32 23
矩阵元素
矩阵元素可以是任何matlab表达 式 ，可以是实数 ，也可以是复 数，复数可用特殊函数I，j 输入 a=[1 2 3;4 5 6] x=[2 pi/2;sqrt(3) 3+5i]
《MATLAB数值计算》PPT课件
ans =
-5.18325528043789
2.17062070347062
-0.83694739215044
0.84958196911772
注意：在上面的程序中，数字格式都设为长(long)型，若改
为短(short)型，结果会有差别，
根据需要可执行 MATLAB 窗口的 Fle | Preferences命令进
第3章 MATLAB数值计算
20.01.2021
精选课件ppt
1
第3章 MATLAB数值计算
3.1 多项式 3.2 插值和拟合 3.3 数值微积分 3.4 线性方程组的数值解 3.5 稀疏矩阵 3.6 常微分方程的数值解
精选课件ppt
2
3.1 多项式
3.1.1 多项式的表达和创建
表示成向量的形式，系数按降序排列 例如
精选课件ppt
11
3.2 插值和拟合
3.2.1 多项式插值和拟合 ➢插值
已知 节点
构造函数
使得
精选课件ppt
12
➢拟合
拟合就是要找出一个曲线方程式（多项式拟合就是设 法找一个多项式），使得它与观测数据最为接近，这时 不要求拟合多项式通过全部已知的观测节点。
1．多项式插值函数(interp1)
yi = interp1(x,y,xi,method) 对应于插值函数
31
精选课件ppt
6
【例 3.4】 利用 polyval找出多项式 在[-1,4]间均匀分布的 5个离散点的值。 >> x=linspace(-1,4,5) % 在[-1,4]区间产生5个离散点
>> p=[1 4 7 -8]; >> v=polyval(p,x) x=
matlab数值计算
第三章数值分析一．多项式例ex3_1.m 多项式的定义、求根、求导%多项式的显示、系数矩阵及多项式的值disp('y(x)=x^3+5*x^2-9*x+3') %显示多项式表达式p=[1 5 -9 3]; %多项式系数矩阵x=1; %x赋初值y1=polyval(p,x); %计算x点处多项式的值%多项式求根、求导r1=roots(p); %求多项式的根p1=poly(r1); %用根构造多项式dy1=polyder(p); %对多项式求导数%显示结果disp(p); %显示多项式系数矩阵disp(y1); %显示x点处多项式的值disp(r1); %显示多项式的根disp(p1); %显示用根构造多项式的系数矩阵disp(dy1); %显示多项式一阶导数例ex3_2.m 多项式的乘、除disp('a(x)=x^3+2*x^2+3*x+4');disp('b(x)=x^3+4*x^2+9*x+16');a=[1 2 3 4]; %多项式a系数矩阵b=[1 4 9 16]; %多项式b系数矩阵c=conv(a,b) %两个多项式相乘，实际是求两个向量的卷积(Convolution):%∑--+=kiikbiakc1)1()()([d,r]=deconv(c,b) %多项式的除法，实际是卷积的逆运算%(Deconvolution), d,r为商多项式与余多项式:%∑=-+= -kiikqjakrkc1)1()()()(例ex3_3.m 多项式拟合x=0:0.1:1;y=[2.1,2.3,2.5,2.9,3.2,3.3,3.8,4.1,4.9,5.4,5.8];n=5; %拟合多项式的阶数取5p=polyfit(x,y,n); %用5阶多项式拟合x、y向量给定的数据y1=polyval(p,x); %计算x点处多项式的值plot(x,y,'o',x,y1,'-'); %绘x、y点图和拟合后的x、y1点图legend('y(x)','y1(x)'); %图例标注例ex3_4.m阶数对拟合效果的影响x=0:0.5:10;y=sqrt(x)+3*sin(x);n=2;p=polyfit(x,y,n);p2=polyval(p,x);n=4;p=polyfit(x,y,n);p4=polyval(p,x);n=6;p=polyfit(x,y,n);p6=polyval(p,x);n=8;p=polyfit(x,y,n);p8=polyval(p,x);plot(x,y,'o',x,p2,'-',x,p4,'-',x,p6,'-',x,p8,'-'); legend('y(x)','n=2','n=4','n=6','n=8');二、插值例ex3_5.m一维插值x=0:1:2*pi;y=sin(x);xi=0:0.1:6.5;yi1=interp1(x,y,xi,'nearst'); %最临近插值yi2=interp1(x,y,xi,'linear'); %线形插值yi3=interp1(x,y,xi,'cubic'); %三次多项插值yi4=interp1(x,y,xi,'spline'); %三次样条插值plot(x,y,'o',xi,yi1,xi,yi2,xi,yi3,xi,yi4);legend('y=sin(x)','nearst','linear','cubic','spline'); 例ex3_6.m高维函数插值（不讲！）[x,y]=meshgrid(-3:0.3:3);z=peaks(x,y);[xi,yi]=meshgrid(-3:0.1:3);zi=interp2(x,y,z,xi,yi,'spline');%zi=interp2(x,y,z,xi,yi,'cubic');%zi=interp2(x,y,z,xi,yi,'linear');%zi=interp2(x,y,z,xi,yi,'nearest');%surf(x,y,z);mesh(x,y,z);hold on;%surf(xi,yi,zi+15);mesh(xi,yi,zi+15);hold off;axis tight;三、快速富叶变换与逆变换（不讲！）例ex3_7.m快速富叶变换t=0:0.001:0.6;x=sin(2*pi*50*t)+sin(2*pi*120*t);y=x+2*randn(size(t));subplot(2,1,1);plot(y(1:50));%---------------y=fft(y); %快速富叶变换f=0:500;subplot(2,1,2);plot(f,y(f+1));例ex3_8.m滤波clear;x=linspace(0,2*pi,64);s=5*sin(x)+2*sin(5*x)+randn(size(x));f=fft(s);f1(1:9)=f(1:9);f1(56:64)=f(56:64);s1=ifft(f1);subplot(2,2,1);plot(x,s);%---------------subplot(2,2,3);fx=0:63;plot(fx,f(fx+1));%---------------subplot(2,2,2);plot(fx,f1(fx+1));%---------------subplot(2,2,4);plot(x,s1);四、稀疏矩阵（不讲！）例ex3_9.ma=[0 0 0 50 2 0 01 3 0 00 0 4 0];s=sparse(a) %转为稀疏矩阵形式b=full(s) %转为全元素矩阵形式c=sparse([3 2 3 4 1],[1 2 2 3 4],[1 2 3 4 5],4,4)%创建稀疏矩阵，c=[I,J,S,m,n,] I、J--行下标、列下标,S—按列排列的所有非零元素构成%的向量,m、n—待生成的稀疏矩阵行、列维数例ex3_10.mc=sparse([3 2 3 4 1],[1 2 2 3 4],[1 2 3 4 5],4,4)nnz(c) %稀疏矩阵的非零元素总数nonzeros(c) %稀疏矩阵的非零元素数值[i,j,s]=find(c) %找出稀疏矩阵的所有非零元素，按列排列[n,m]=size(c) %稀疏矩阵行、列维数d=c+ones(4,4) %稀疏矩阵加1矩阵五、数值积分例ex3_11.m(1)建立函数fn1function y=fn1(x)y=exp(-x.*x)(2)对被积函数fn1进行积分s1=quad('fn1',0,1,0.001); %采用Simpson法计算积分%0,1: 下限、上限%0.001: 精度s2=quad8('fn1',0,1,0.001,1);%采用8样条Newton-Cutes%公式求数值积分%最后的1:显示积分过程，0：不显示积分过程六、微分方程的数值解先将高阶微分方程降阶处理，转换为一阶微分方程组，再利用ode45函求解。
matlab第3章
第7章MATLAB科学计算¾方程求解¾概率统计¾插值、拟合¾数值微积分¾最优化求解其它常用的matlab 数值计算命令¾max,min¾mean, median¾sum 求和, prod 求积¾cumsum 求和, cumprod 求积¾std 标准方差¾corrcoef 相关系数¾sort 元素排序¾离散傅里叶变换fft,fft2,fftn__ifft第7章MATLAB 数值计算作业¾1.编写傅立叶变换的matlab 程序与matlab 自带的fft 进行比较，并分析冲击信号的傅立叶变换。
（若不了解冲击信号，可计算方波的傅里叶变换，方波幅度为1，周期为10，方波个数为10，占空比为0.5）。
∑=−−−=Nm Nk m j em f k F 1/)1)(1(2)()(π编写的DFT 函数：function X=mydft(x )N = length(x );W=exp(-2*i*pi/N);X=zeros(1,N);for k=1:NX(k )=sum(x .*W.^((0:N -1)*(k -1)));end∑=−−−=N m N k m j em f k F 1/)1)(1(2)()(πx = [0 0 0 0 0 1 1 1 1 1]; X = [x x x x x x x x x x]; y = mydft(X);plot(abs(y))y1=fft(X);plot(abs(y1));¾y = fftshift(mydft(X));¾>> plot(abs(y))第3章MATLAB符号计算¾Maple优势在于符号运算，¾Mathematic符号运算和数值计算均不差，图像处理或者数据可视化较差¾Matlab强项是数值计算和数据可视化，¾MathCAD各方面均弱一些，但易学。
matlab第三章MATLAB数值计算
min函数的用法和max完全相同。
例 分析下列程序的功能
x=[4 5 6;1 4 8]; y=[1 7 5;4 5 7]; p=max(x,y) ; p
分析：取两个2×3的二维数组x和y同一位置上的元素值
大者构成一个新矩阵p。
2 平均值和中值
求数据序列平均值的函数是mean，求数据序列中值
的函数是median。
4
累加和与累乘积
向量的累加和与累乘积
cumsum(X)
返回向量X累加和向量。 cumprod(X) 返回向量X累乘积向量。 矩阵的累加和与累乘积
cumsum(A)
返回一个矩阵，其第i列是A的第i列的累加和向量。
cumprod(A) 返回一个矩阵，其第i列是A的第i列的累乘积向量。 cumsum(A,dim) 当dim为1时，该函数等同于cumsum(A)；当dim为2时，
(2) 对角阵
X = diag(v,k) 当v是n个元素的向量时，返回有第k个对角线的 n+abs(k)顺序的方阵，k = 0（可省略）代表主对角线， k > 0代表上方的次对角线， k < 0代表下方的次对角
y3=median(x,2); %按行方向,求矩阵的中值
y4=mean(x); y5=mean(x,1); %按列方向,求矩阵的平均值
y6=mean(x,2); %按行方向,求矩阵的平均值
3 求和与求积
向量的求和与求积
sum(X) 返回向量X各元素的和。 prod(X) 返回向量X各元素的乘积。
3.3 矩阵操作 3.3.1 矩阵的结构变换
(1) 转置
转置运算的操作符： ′ 求A的转置，运算表达式为 A′ 其中A可以是行向量、列向量和矩阵。 若矩阵中的元素为复数， 则′对矩阵元素做转置共轭； 若仅对矩阵中元素做转置，则 可用 .′操作符
实验3 MATLAB 数值计算(2)
实验3 MATLAB 数值计算（2）
目的和要求：
（1）了解多项式的运算。
（2）熟练掌握MATLAB二维曲线的绘制。
内容和步骤：
1．多项式的运算式的运算
（1）多项式的运算。
已知表达式G(x)=(x-4)(x+5)(x2-6x+9), 展开多项式形式；求导；并计算当x在[0，20]范围变化时G(x)的值；计算出G(x)=0的根。
多项式相乘conv;
求导polyder;
计算零点，即求根roots
解：
展开为多项式
求导：
求零点：
（2）多项式的拟合和插值。
x在[0，20]范围内，计算多项式y=x4-5x3-17x2+129x-180 的值y；并根据x和y进行二阶、三阶和四阶拟合。
并绘出拟合曲线。
对多项式y进行插值，计算在5.5处的值。
多项式拟合p=ployfit(x,y,n)
插值yi=interp1(x,y,xi,’method’)
2．绘制二维曲线
绘制的图形窗口分割为一行两列，窗口左上角画一正弦曲线，y=sin(2t),t:[0.2π];窗口右上角画3条衰减的单边指数曲线y=e-t, y=e-2t,和y=e-3t, t:[0,2]。
在图上添加标题，将3条曲线用不同的线型，并添加图例。
MATLAB数值运算.pdf
第3章 MATLAB 数值运算教学提示：每当难以对一个函数进行积分或者微分以确定一些特殊的值时，可以借助计算机在数值上近似所需的结果，从而生成其他方法无法求解的问题的近似解。
这在计算机科学和数学领域，称为数值分析。
本章涉及的数值分析的主要内容有插值与多项式拟合、数值微积分、线性方程组的数值求解、微分方程的求解等，掌握这些主要内容及相应的基本算法有助于分析、理解、改进甚至构造新的数值算法。
教学要求：本章主要是让学生掌握数值分析中多项式插值和拟合、牛顿-科茨系列数值求积公式、3种迭代方法求解线性方程组、解常微分方程的欧拉法和龙格-库塔法等具体的数值算法，并要求这些数值算法能在MATLAB 中实现。
3.1 多 项 式在工程及科学分析上，多项式常被用来模拟一个物理现象的解析函数。
之所以采用多项式，是因为它很容易计算，多项式运算是数学中最基本的运算之一。
在高等数学中，多项式一般可表示为以下形式：120121()n n n n n f x a x a x a x a x a −−−=+++++…。
当x 是矩阵形式时，代表矩阵多项式，矩阵多项式是矩阵分析的一个重要组成部分,也是控制论和系统工程的一个重要工具。
3.1.1 多项式的表达和创建在MATLAB 中，多项式表示成向量的形式，它的系数是按降序排列的。
只需将按降幂次序的多项式的每个系数填入向量中，就可以在MATLAB 中建立一个多项式。
例如，多项式43231529s s s s +−−+在MATLAB 中，按下面方式组成一个向量x = [1 3 -15 -2 9]MATLAB 会将长度为n +1的向量解释成一个n 阶多项式。
因此，若多项式某些项系数为零，则必须在向量中相应位置补零。
例如多项式41s +在MATLAB 环境下表示为y = [1 0 0 0 1]3.1.2 多项式的四则运算多项式的四则运算包括多项式的加、减、乘、除运算。
下面以对两个同阶次多项式MATLAB 基础及其应用教程·66··66·32()234a x x x x =+++，32()4916b x x x x =+++做加减乘除运算为例，说明多项式的四则运算过程。
MATLAB 数值计算(2)
数值计算MATLAB 数值计算第四章MATLAB1主要内容基本数据运算数据统计与分析数据插值与曲线拟合多项式计算数值微积分线性方程组求解非线方与问求解非线性方程与最优化问题求解2常微分方程的数值求解多项式计算N次多项式表示为–P(x)=a 0x n +a 1x n-1+a 2x n-2…a n-1x+a nMatlab中n次多项式用一个长度为n+1的行向量（系数向量）表示–[a a …a n-1a [01n 1n ]多项式的四则运算–多项式的加减运算»系数向量的加减运算要求次数相同不足时用“»要求次数相同，不足时用“0”补起——向量化处理»例54322()352756f x x x x x x =−+−++3()353g x x x =+−多项式乘法运算–函数conv(P1,P2)用于求多项式P1和P2的乘积。
这里，P1、P2是两个多项式系数向量 多项式除法运算–函数[Q,r]=deconv(P1,P2)用于对多项式P1和P2作除法运算。
其中Q 返回多项式P1除的商式的余式这里仍是多项式系数向量以P2的商式，r 返回P1除以P2的余式。
这里，Q 和r 仍是多项式系数向量。
–deconv 是conv 的逆函数，即有P1=conv(P2,Q)+r 。
5432−− 例–求f(x)+g(x)、f(x)-g(x)。
2()352756()353f x x x x x xg x x x =+++=+−–求f(x)×g(x)、f(x)/g(x)。
–f=[3,-5,2,-7,5,6];g=[3,5,-3];g1=[0,0,0,g];–f+g1%求f(x)+g(x)f+g1 %求f(x)+g(x)–f-g1 %求f(x)-g(x)–conv(f,g) %求f(x)*g(x)[]()求()/()商式送余式送4–[Q,r]=deconv(f,g) %求f(x)/g(x)，商式送Q，余式送r。
《现代数值计算方法(MATLAB版)》习题解答
= 0 0 0 0 0
2.7 提示: Bs = (D − L)−1 U = − 1 2 0 值 λ1 = 0, λ2 = λ3 0 1 BJ = 2 −2 1 Jacobi 迭代发散. = −1 , 2
1 2
0 1
1 2
2.2218 ≤ n ≤ 2.9208 ⇒ n √ = 2. 1.8 提示: x1,2 =
282 − √781 28+ 783
= 28 ±
√
783, x1 = 28 + 27.982 = 55.982 ≈ 55.98, x2 = 28 −
1−cos2 1◦ 1+cos 1◦
=
1 55.982
≈ 0.01786. =
√
5 2
> 1, 故
2.8 提示: (1) A = 1 3 a > 1, ⇒ a3 − 14a + 12 > 0, Seidel 迭代收敛.
a > 0, a 2 − 1 > 0, ⇒ 2 , 当 |a| > 5 时, Jacobi 迭代收敛. (2) a3 − 14a + 12 > 0, a 所以, 当 a ≥ √ 14 时, A 对称正定, 从而 Gauss-
10 +1+10
1.11 (1) (A) 比较准确; (2) (A) 比较准确. 1.12 算法 2 准确. 在算法 1 中, ε0 ≈ 0.2231 带有误差 0.5 × 10−4 , 而这个误差在以后的每次计算中 顺次以 41 , 42 , · · · 传播到 In 中. 而算法 2 中的误差是按
现代数值计算方法(MATLAB版)第3章(1)
akk ,
(k )
(3.5)
福建师范大学
数计学院
8/15
k = 1, · · · , n − 1, mik = aik /akk , aij
(k +1) (k ) (k ) (k )
aik
(k +1) (k )
= 0, bi
(k +1)
= aij − mik akj ,
= bi − mik bk .
Back Close
(k )
(k )
(i = k + 1, · · · , n; j = k + 1, · · · , n.) 3
xn = bn /ann , k = n − 1, · · · , 1, xk =
(k ) bk n (k ) (k )
福建师范大学
(n)
(n)
−
j =k +1
akj xj
··· ··· ···
(1) a1n (2) a2n (2) a3n
(1) b1 (2) b2 (2) b3
福建师范大学
数计学院
6/15
··· ··· ··· ··· an2 an3 · · · ann bn
(2) (2) (2) (2)
Back Close
(2) aij
(1)
,
A(k) A(k) k
Dk = det Dk = 0
A11
(k −1)
A12
(k −1) (k )
(k −1) (k ) = a(1) 11 · · · ak −1,k −1 akk .
akk
akk = 0.
(k )
MATLAB第三章数值数组及其运算
行向量
如：array=[2, pi/2, sqrt(3), 3+5i]
x=[1,2,3,4,5都已知.如对 少量实验数据的处理可用此种方法.
4
(2) 冒号生成法： array=a: inc: b
<向量>
a---数组的第一个元素
inc---采样点之间的间隔, 即步长. 最后一个元素不一定等于b, 其大小为b’=a +inc*[(b-a)/inc]; 步长可以省略, 默认为 1; inc可以取正数或负数, 但要注意当取正时,要保证b>a, 数 组最后一个元素不超过b, 取负时b<a, 最后一个元素不小于b.
(2) 数值计算解法
delt=0.01; x=0:delt:4;
y=exp(-sin(x));
sx=delt*cumtrapz(y);
plot(x,y, 'r', 'LineWidth', 6); hold on;
plot(x, sx, '.b', 'MarkerSize', 15);
plot(x, ones(size(x)), 'k');
a inc>0 b
b inc>0 a
特点: 等差数列
方便对数据之间的间隔(步长)进行控制.但要注意三个数值之 间的关系，可能得到空数组．另外要注意生成的数组的元素的 个数．如x=a: (b-a)/n :b (b>a)得到n+1个元素的数组.
5
x=1:5x=[1,2,3,4,5]
y=5:-1:1y=[5, 4, 3, 2, 1]
2. 在命令窗中输入MyMatrix
11
3.5 二维数组的标识 (mxn, m>1, n>1)
数值计算与MATLAB方法课后答案
第一章习题1. 序列满足递推关系，取及试分别计算，从而说明递推公式对于计算是不稳定的。
n1 1 0.01 0.00012 0.01 0.0001 0.0000013 0.0001 0.000001 0.000000014 0.000001 0.0000000110-105 0.00000001 10-10n1 1.000001 0.01 0.0000992 0.01 0.000099 -0.000099013 0.000099 -0.00009901-0.010000994 -0.00009901 -0.01000099-1.00015 -0.01000099-1.0001初始相差不大，而却相差那么远，计算是不稳定的。
2. 取y0=28，按递推公式，去计算y100，若取（五位有效数字），试问计算y100将有多大误差？y100中尚留有几位有效数字？解：每递推一次有误差因此，尚留有二位有效数字。
3．函数，求f(30)的值。
若开方用六位函数表，问求对数时误差有多大？若改用另一等价公式计算，求对数时误差有多大？设z=ln(30-y)，，y*， |E(y)| 10-4z*=ln(30-y*)=ln(0.0167)=-4.09235若改用等价公式设z=-ln(30+y)，，y*， |E(y)|⨯10-4z*=-ln(30+y*)=-ln(59.9833)=-4.094074．下列各数都按有效数字给出，试估计f的绝对误差限和相对误差限。
1)f=sin[(3.14)(2.685)]设f=sin xyx*=3.14, E(x)⨯10-2, y*=2.685, E(y)⨯10-3,sin(x*y*)=0.838147484, cos(x*y*)=-0.545443667⨯(-0.5454) ⨯⨯10-2+3.14(-0.5454) ⨯⨯10-3|⨯10-2⨯10-2|E r(f)| ⨯10-2⨯10-2<10-22)f=(1.56)设f = x y ,x*=1.56, E(x)⨯10-2, y*=3.414, E(y)⨯10-3,⨯⨯⨯10-2⨯⨯⨯10-3|⨯⨯⨯10-2⨯⨯⨯10-3|=0.051|E r(f)| =0.01125．计算，利用下列等式计算，哪一个得到的结果最好，为什么？6．下列各式怎样计算才能减少误差？7. 求方程x2-56x+1=0的二个根，问要使它们具有四位有效数字，至少要取几位有效数字？如果利用伟达定理， 又该取几位有效数字呢？解一:若要取到四位有效数字,如果利用伟达定理,解二:由定理二,欲使x1,x2有四位有效数字,必须使由定理一知，∆至少要取7位有效数字。
MATLAB 数学实验 第三章
微积分符号计算 diff(f) — 对缺省变量求导数 diff(f,v) — 对指定变量 v 求导数 diff(f,v,n) —对指定变量 v 求n阶导数 对指定变量 阶导数 int(f) — 对f表达式的缺省变量求积分 表达式的缺省变量求积分 int(f,v) — 对f表达式的 变量求积分 表达式的v变量求积分 表达式的 int(f,v,a,b) — 对f表达式的 变量在 b] 表达式的v变量在 表达式的 变量在[a, 区间求定积分
绕X轴旋转的旋转曲面体积 轴旋转的旋转曲面体积 2π V = π ∫ [ f ( x )]2 dx 0 syms a b x f=exp(a*x)*sin(b*x); f1=subs(f,a,-0.2); f2=subs(f1,b,0.5); V=pi*int(f2*f2,x,0,2*pi) double(V) V =pi*(-125/116*exp(-4/5*pi)+125/116) ans = 3.1111
16/20
1 2 y2 例3.26 解微分方程 y ′ = 2 x +1
y ( 0) = 0
命令格式:dsolve('eq1',,'con1',,'x') y的一阶导数—— Dy, y的二阶导数—— D2y
y = dsolve('Dy=1/(1+x^2)-2*y^2','y(0) = 0','x') y= 2*x/(2*x^2+2) 符号解: 符号解: y(x)= x / (1 + x 2)
12/20
定积分数值计算命令 quad(f, a, b) t 例3.14 计算积分上限函数值 F (t ) = ∫
MATLAB数值计算
自然对数〔以e为底〕
常用对数〔以10为底〕
12
2.2 矩阵和数组
2.2.1 矩阵的赋值 2.2.2 向量的生成 2.2.3 矩阵元素 2.2.4 复数表示
13
2.2.1 矩阵的赋值
(1). 直接输入法创立矩阵
矩阵的全部元素必需放在方括号“[]”内；
矩阵元素之间必需用逗号“，”或空格隔开；
30
〔2〕.矩阵的赋值
全下标方式：A(i,j)=B给A矩阵的局部元素赋值 则B矩阵的行列数必需等于A矩阵的行列数。
A(1:2,1:3)=[1 1 1;1 1 1] A=1 1 1 6 2
1 1147 7 5715 0 3454 23 13 6 0 3
31
单下标方式：A(s)=b,b为向量，元素个数必需等 于A矩阵的元素个数。
44
〔6〕. 矩阵的翻转
矩阵的翻转及对角化操作函数
命令
说明
flipud(A) 矩阵作上下翻转
fliplr(A) rot90(A)
矩阵作左右翻转 矩阵逆时针翻转90°
diag(A) 提取矩阵A的对角元素，返回列向量
diag(V) 以列向量V作对角元素创建对角矩阵
tril(A)
提取矩阵A的下三角矩阵
triu(A)
61
关系运算的规章
当参与运算的矩阵是两同维矩阵A和B时， 关系运算的结果是将矩阵A 和B 下标一样的对 应元素逐一进展关系比较，假设关系成立则比 较结果值为“1”，假设关系不成立则比较结 果值为“0”。也即关系运算的结果是生成一 个与A 和B 维数一样的矩阵，其元素值为“0” 或“1”。 算术运算比关系运算具有更高的优先权。
取值 用于结果的缺省变量名 圆周率 计算机的最小数，当和1相加就产生一个比1大的数 浮点运算数 无穷大，如1/0 不定量，如0/0
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数计学院
16/16
ak ck−1 |¯ bk | = bk − ck−1 ≥ |bk | − |ak | · > |bk | − |ak | > 0, bk−1 bk−1 bk = 0, k = 2, · · · , n. , .
Back Close
福建师范大学
数计学院
10/16
MATLAB >> A=[2 -1 4 -3 1;-1 1 2 1 3;4 2 3 3 -1; -3 1 3 2 4; 1 3 -1 4 4]; >> b=[11 14 4 16 18]’; >> x=magauss2(A,b); x’
Back Close
x = 1.0000 §3.3 , b1 c1 a2 b 2 c 2 ... ... ... an−1 bn−1 cn−1 an b n Gauss x1 x2 . . . xn−1 xn d1 d2 = . . . dn−1 dn . 2.0000 1.0000 -1.0000 4.0000
(3.8)
福建师范大学
数计学院
13/16
(3.9)
Back Close
function x= machase(a,b,c,d) % % % % : : x= machase(a,b,c,d) b d ,x , c Ax=d a ,
14/16
,
福建师范大学
数计学院
n=length(a); for k=2:n b(k)=b(k)-a(k)/b(k-1)*c(k-1); d(k)=d(k)-a(k)/b(k-1)*d(k-1); end x(n)=d(n)/b(n); for k=n-1:-1:1 x(k)=(d(k)-c(k)*x(k+1))/b(k);
Back Close 9/16
福建师范大学
数计学院
3.5 2 −1 4 −3 1 −1 1 2 1 3
magauss2.m 4 2 3 3 −1 −3 1 3 2 4 x1 1 3 x2 −1 x3 4 x4 x5 4 11 14 = 4 16 18
3/16
−10000.0 −10000.0
x = 1.0, 0.0001x + 1.0x = 1.0, 1 2 2 =⇒ x1 = 0.0. −10000.0x2 = −10000.0.
Back Close
, , , , ,
. 0.0001 10000 . .
(3.6)
福建师范大学
rk
k
.
数计学院
6/16
Gauss A,
) b, k := 1
k = 1, · · · , n − 1 , rk ,
(k ) k ≤i≤n
ark k = max |aik |, = 0, , (A(k), b(k)) , k, rk
Back Close
(2) (3)
rk > k ,
Back Close 8/16
% m=A(k+1:n,k)/A(k,k); A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m*A(k,k+1:n); b(k+1:n)=b(k+1:n)-m*b(k); A(k+1:n,k)=zeros(n-k,1); if flag~=0, Ab=[A,b], end end % x=zeros(n,1); x(n)=b(n)/A(n,n); for k=n-1:-1:1 x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k); end
Back Close
3.7
(3.7)
:
|bi| > |ci| > 0, |bi+1| > |ai+1| > 0, i = 1, · · · , n − 1,
福建师范大学
. (3.8)-(3.9) ¯ b1 = b1 = 0. k≥2 , , ¯ bk = 0, (k = 1, 2, · · · , n) .
r 1 ↔r 2 − − − →
1.0
1.0 2.0
Back Close
, 1.0 − 0.0002 = 1.0.
1.0 − 0.0001 = (0.1 − 0.00001)101 (
),
福建师范大学
数计学院
x = 1.0, 1.0x + 1.0x = 2.0, 2 1 2 =⇒ x1 = 1.0. 1.0x2 = 1.0. , . A(1) = A, ar 1 1 , ar11 = max |ai1 |
Back Close
end 3.6 machase.m
福建师范大学
数计学院
4 1 x1 1 4 1 x2 1 4 1 x3 . . . . . . . 1 . . ... 4 1 x 49 1 4 x50 MATLAB
5 6 6 = . . . 6 5
15/16
>> a=ones(50,1); b=4*ones(50,1); c=ones(50,1); >> d=6*ones(50,1); d(1)=5; d(50)=5; >> x=machase(a,b,c,d) x = (x1, x2, · · · , x50) = (1.0, 1.0, · · · , 1.0).
(k )
2/16
akk = 0,
(k )
x1 =
x2 =
Back Close
4 4 0.1)105 = (0.0000 − 0.1)105 = −10000.0 (
. 1.0 − 10000.0 = (0.00001 − ), , 2.0 − 10000.0 =
福建师范大学
数计学院
−10000.0, 0.0001 1.0 1.0 0.0001 1.0 1.0 r2 −104 r1 − − − − − → 1.0 1.0 2.0 0 1.0 − 10000.0 2.0 − 10000.0 − − − → 0.0001 0 1.0 1.0 .
d1 d2 . . . dn−1 dn ¯1 d ¯2 d . . . ¯n−1 d ¯n d
福建师范大学
数计学院
12/16
Back Close
¯ ¯1 = d1, b1 = b1 , d ak ¯ ck−1, (k = 2, · · · , n) bk = bk − b k −1 ¯k = dk − ak d ¯k−1 d bk−1 ¯n d xn = ¯ , bn ¯k − ck xk+1 d xk = , k = n − 1, · · · , 2, 1. ¯ bk , , MATLAB • MATLAB %machase.m 6n − 5 . ,
福建师范大学
数计学院
11/16
(3.7)
Back Close
b c 1 1 a2 b2 c2 ... ... ... an−1 bn−1 cn−1 an b n ¯ b c 1 1 ¯ b2 c 2 ... ... → ¯ bn−1 cn−1 ¯ bn
aik
(k +1) (k )
= 0, bi
(k +1)
= aij − mik akj ,
= bi − mik bk .
福建师范大学
(k )
(k )
数计学院
7/16
−
j =k +1
akj xj
(k )
akk .
(k )
MATLAB MATLAB
function x=magauss2(A,b,flag) % : Gauss Ax=b
4/16
,
福建师范大学
数计学院
0.0001 1.0 1.0 1.0
1.0 2.0 0.0001 1.0 1.0 1.0 1.0 2.0 r2 −0.0001r1 − − − − − − → 0 1.0 − 0.0001 1.0 − 0.0002 1.0 1.0 2.0 . − − − → 0 1.0 1.0
1≤i≤n (1) (1) (1)
5/16
. . Gauss 1 , Gauss 1
.
r1 > 1,
r1
1
.
Back Close
,
Gauss
(k )
k
k ≤i≤n
,
(k )
ark k = max |aik | . rk > k , Gauss 3.2 ( 1 2 (1)
(k ) (k ) ar k k
Back Close
% % %
ห้องสมุดไป่ตู้
: x=magauss(A,b,flag), A flag=0, 0, x ,
, b
, ,
福建师范大学
数计学院
if nargin<3,flag=0;end n=length(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; t=b(k); b(k)=b(p); b(p)=t; end
福建师范大学