现代数值计算方法(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可以是行向量、列向量和矩阵。 若矩阵中的元素为复数， 则′对矩阵元素做转置共轭； 若仅对矩阵中元素做转置，则 可用 .′操作符
- 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
福建师范大学