数值计算
MATLAB 数值计算
第四章MATLAB
1
主要内容
基本数据运算 数据统计与分析 数据插值与曲线拟合 多项式计算 数值微积分 线性方程组求解
非线方与问求解 非线性方程与最优化问题求解2
常微分方程的数值求解
多项式计算
N次多项式表示为
–P(x)=a 0x n +a 1x n-1+a 2x n-2…a n-1x+a n
Matlab中n次多项式用一个长度为n+1的行向量（系数向量
）表示
–[a a …a n-1a [01n 1n ]
多项式的四则运算
–多项式的加减运算
?系数向量的加减运算
要求次数相同不足时用“?要求次数相同，不足时用“0”补起——向量化处理?例
54322()352756f x x x x x x =?+?++3
()353
g 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()353
f x x x x x x
g 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。
多项式的导函数
–p=polyder(P)：求多项式P的导函数
l d(P)求多项式
多项式求值
–MATLAB 提供了两种求多项式值的函数：polyval 与polyvalm ，它们的输入参数均为多项式系数向量P 和自变量x 。–两者的区别在于前者是代数多项式求值，而后者是矩阵多项式求值。
代数多项式求值
–Y=polyval(P,x)
–若为一数值，则求多项式在该点的值；若x 为数值，则求多项式在该点的值；若x 为向量或矩阵，则对向量或矩阵中的每个元素求其多项式的值。
矩阵多项式求值
–polyvalm 函数用来求矩阵多项式的值，其调用格式与polyval 相同，但含义不同。–polyvalm 函数要求x 为方阵，它以方阵为自变量求多项式的值。
polyval 与polyvalm 的区别
–设A 为方阵，P 代表多项式x 3-5x 2+8，那么polyvalm(P,A)的含义是：A*A*A-5*A*A+8*eye(size(A))
6
–polyval(P,A)的含义是：A.*A.*A-5*A.*A+8*ones(size(A))
以多项式x4+8x3-10为例，取一个2×2矩阵为自变量分别用polyval和polyvalm计算该多项式的值。
l l l l计算该多项式的值
–A=[1,8,0,0,-10]; % 多项式系数
–x=[-1,1.2;2,-1.8] % 给出一个矩阵x
–y1=polyval(A,x) % 计算代数多项式的值
–y2=polyvalm(A,x) % 计算矩阵多项式的值
7
多项式求根
n次多项式具有n个根，当然这些根可能是实根，也可能含有若干对共轭复根
–MATLAB提供的roots函数用于求多项式的全部根，其调用格式为：x=roots(P)
?其中P为多项式的系数向量，求得的根赋给向量x，即x(1),x(2),…,x(n)分别代表多项式的n个根。
–P=poly(x)
?若x为具有n个元素的向量，则poly(x)建立以x为其根的多项式，且将该多项式的系数赋给向量P。
8
例
已知5
2.7543)(235+??+=x x x x x f (1) 计算f(x)=0 的全部根。
(2)f(x)=0的根构造一个多项式(2) 由方程f(x)0的根构造个多项式g(x)，并与f(x)进行对比。命令如下命令如下：P=[3,0,4,-5,-7.2,5];
X=roots(P) %求方程f(x)=0的根G=poly(X)%G=poly(X) %求多项式g(x)
9
主要内容
基本数据运算 数据统计与分析 数据插值与曲线拟合 多项式计算 数值微积分 线性方程组求解
非线方与问求解 非线性方程与最优化问题求解10
常微分方程的数值求解
数值微分的两种方法
g()对()进行逼近插值或拟合，然后用逼近函数g()–用多项式或样条函数g(x)f(x)进行逼近（插值或拟合），然后用逼近函数g(x)在点x处的导数作为f(x)在点x处的导数
–用f(x)在点x处的某种差商作为其导数
在MATLAB中，没有直接提供求数值导数的函数，只有计算向前差分的函数diff，其调用格式为：
–DX=diff(X)：计算向量X的向前差分，DX(i)=X(i+1)-X(i)，i=1,2,…,n-1
–DX=diff(X,n)：计算X的n阶向前差分。例如，diff(X,2)=diff(diff(X))
–DX=diff(A,n,dim)：计算矩阵A的n阶差分，dim=1时(缺省状态)，按列计算差分；dim=2，按行计算差分。
例设x由[0,2π]间均匀分布的10个点组成，求sin x的1~3阶差分–X=linspace(0,2*pi,10);
–Y=sin(X);
–DY=diff(Y); %计算Y的一阶差分
–D2Y=diff(Y,2); %计算Y的二阶差分，也可用命令diff(DY)计算
–D3Y=diff(Y,3); %计算Y的三阶差分，也可用diff(D2Y)或diff(DY,2)
11
数值积分基本原理
–求解定积分的数值方法多种多样，如简单的梯形法、辛普生(Simpson)法、牛顿－柯(p )特斯(Newton-Cotes)法等都是经常采用的方法。它们的基本思想都是将整个积分区间[a,b]分成n 个子区间[x i ,x i+1]，i=1,2,…,n ，其中x 1=a ，x n+1=b 。这样求定积分问题就分解为求和问题。
数值积分的实现
–变步长辛普生法
[I,n]=quad(‘fname ’,a,b,tol,trace)
其中fname是被积函数名。a和b分别是定积分的下限和上限。tol用来控制积分精度，缺省时取tol=0.001。trace控制是否展现积分过程，若取非0则展现积分过程，取0则不展现，缺省时取trace=0。返回参数I即定积分值，n为被积函数的调用次数牛顿柯斯特法
–牛顿-柯斯特法[I,n]=quad8(‘fname ’,a,b,tol,trace)–Lobbato法——高精度[I n]=quadl(‘’a b tol trace)
[I,n]quadl(fname ,a,b,tol,trace)–被积函数由一个表格定义
trapz(X,Y) 其中向量X 、Y 分别代表样本点和样本值，定义函数关系Y=f(X)二重定积分
12
–I=dblquad(’fname ’,a,b,c,d,tol,trace)
fname为描述被积函数的字符串变量，定义有两种方式–建立关于被积函数的函数文件
function ex=ex(x)
ex=exp(-x.^2);
e e p(^2);
–使用语句函数(内联函数) Inline()
?inline()，第一个变量为被积函数本身，第二个变量为自变量，也可
i li()第个变量为被积函数本身第二个变量为自变量也可
以带多个自变量
g=inline('exp(-x.^2)'); %定义一个语句函数g(x)=exp(-x^2)
13
例
用两种不同的方法求定积分–方法一：先建立一个函数文件ex.m ：∫?=
1
2
dx
e
I x 方法先建个函数文件function ex=ex(x)ex=exp(-x.^2);
然后在命令窗口输入命令：MATLAB 命令窗口，输入命令：format long
I=quad('ex',0,1) %注意函数名应加字符引号I =
0.74682418072642I=quadl('ex',0,1)I =
0.74682413398845
–方法二：也可不建立关于被积函数的函数文件，而使用语句函数(内联函数)求解命令如
，命令如下：g=inline('exp(-x.^2)'); %定义一个语句函数g(x)=exp(-x^2)I=quadl(g,0,1) %注意函数名不加引号14I =
0.74682413398845
format short
主要内容
基本数据运算 数据统计与分析 数据插值与曲线拟合 多项式计算 数值微积分 线性方程组求解
非线方与问求解 非线性方程与最优化问题求解15
常微分方程的数值求解
线性方程组求解
对于线性方程组AX=b ，设A为m×n的矩阵，则根据m和n的关系，可以将方程组分类：
关系可以将方程组分类
–当m = n时，恰定方程组，可以求得精确解
–当m > n时，超定方程组，可以求得最小二乘解。
–当m < n时，欠定方程组，可以求得基础解，该解最多有m个非零元素
16
求解
恰定方程组求解
恰定方程组
方法一：利用左除运算符直接求解：x=A\b
方法二：利用矩阵求逆解法，即x=inv(A)*b
方法三：利用LU分解法求解
–LU分解
分解就是将个矩阵表示为个交换下角矩阵和个角矩?矩阵的LU分解就是将一个矩阵表示为一个交换下三角矩阵和一个上三角矩阵的乘积形式。
?线性代数中已经证明，只要方阵A是非奇异的，LU分解总是可以进行的–MATLAB提供的lu函数用于对矩阵进行LU分解，其调用格式为：
分解其调用格式为
?[L,U]=lu(A)：产生一个上三角阵U和一个变换形式的下三角阵L(行交换)，使之满足A=LU。注意，这里的矩阵A必须是方阵。
?[L,U,P]=lu(A)：产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P，使之满足PA=LU。当然矩阵A同样必须是方阵。
–LU分解后，线性方程组Ax=b的解x=U\(L\b)或x=U\(L\Pb)，这样可以大大提实现分解后线性方程组A b U\(L\b)U\(L\Pb)这样可以大大提
高运算速度。
17
例
用LU 分解求解线性方程组??
?
?=+?=+?+97513
524321x x x x x x x 命令如下：????
=??+=?+0
46624321432421x x x x x x x A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];
b=[13,-9,6,0]';[,,,];[L,U]=lu(A);x=U\(L\b)
18
一般来说，对维数不高，条件数不大的矩阵，前述三种解法所得的结果差别不大。
法所得的结果差别不大
MATLAB中，出于对算法稳定性的考虑，行列式及逆的计算大都在LU分解的基础上进行。
种方法都是在为非奇异矩阵的前提下进行的
三种方法都是在A为非奇异矩阵的前提下进行的。
–当A接近奇异矩阵时，MATLAB将给出警告信息。
–当发现A是奇异矩阵时，计算结果为inf，并且给出警告信息。
当发现A是奇异矩阵时计算结果为inf并且给出警告信息–当矩阵A是病态矩阵时，给出警告信息。
19
例
求解方程组
123123231
4562x x x x x x ++=??
++=?–A = [1, 2, 3; 4, 5, 6; 7, 8, 9];1
237893
x x x ?++=?–b = [1; 2; 3];–x = A\b
–Warning: Matrix is close to singular or badly scaled.
%系统给出警告，结果可能不正确–rank(A)–ans =–
2
–注意：在求解方程时，尽量不要用inv(A)*b命令，而应采用A\b 的解法因为后者的计算速度比前者快精度高尤其当矩阵A 的解法。因为后者的计算速度比前者快、精度高，尤其当矩阵A
的维数比较大时。另外，A\b命令的适用行较强，对于非方阵A
，也能给出最小二乘解。
20
matlab数值计算功能 1,基础运算 （1）多项式的创建与表达 将多项式(x-6)(x-3)(x-8)表示为系数形式 a=[6 3 8] % 写成根矢量 pa=poly(a)% 求出系数矢量 ppa=poly2sym(pa,'x') % 表示成符号形式 ezplot(ppa,[-50,50]) 求3介方阵A的特征多项式 a=[6 2 4;7 5 6;1 3 6 ]; pa=poly(a)% 写出系数矢量 ppa=poly2sym(pa) %表示成符号形式 ezplot(ppa,[-50,50]) % 绘图 求x^3-6x^2-72x-27的根。 a=[1,-6,-72,-85]; % 写出多项式系数矢量 r=roots(a) % 求多项式的根 （2）多项式的乘除运算 c=conv(a,b) %乘法 [q,r]=deconv(c,a)% 除法 求a(s)=s^2+2s+3乘以b(s)=4s^2+5s+6的乘积 a=[1 2 3] b=[4 5 6] % 写出系数矢量 c=conv(a,b) c=poly2sym(c,'s') % 写成符号形式的多项式 展开(s^2+2s+2)(s+4)(s+1)并验证结果被(s+4),(s+3)除后的结果。c=conv([1,2,2],conv([1,4],[1,1])) cs=poly2sym(c,'s') c=[1 7 16 18 8] [q1,r1]=deconv(c,[1,4]) [q2,r2]=deconv(c,[1,3]) cc=conv(q2,[1,3]) test=((c-r2)==cc)
其他常用的多项式运算命令 pa=polyval(p,s) % 按数组规则计算给定s时多项式的值 pm=polyvalm(p,s)% 按矩阵规则计算给定s时多项式的值 [r,p,k]=residue(b,a) % 部分分式展开，b,a分别是分子，分母多项式系数矢量。r,p,k分别是留数，极点和值项矢量。 p=poly(x,y,n) % 用n介多项式拟合给定的数据 polyder(p) %多项式微分 x=[1 2 3 4 5]; y=[5.5 43.1 128 290.7 498.4]; p=polyfit(x,y,3) x2=1:0.1:5; y2=polyval(p,x2); plot(x,y,'o',x2,y2) 2，线性代数 1，求解方程的根 t=[0 0.3 0.8 1.1 1.6 2.3]'; y=[0.82 0.72 0.63 0.60 0.55 0.50]'; e=[ones(size(t)) exp(-t)] c=e\y t1=[0:0.1:2.5]'; y1=[ones(size(t1)),exp(-t1)]*c; plot(t1,y1,'b',t,y,'ro') 2，逆矩阵及行列式 inv(a) det(a) pinv(a) 3，矩阵分解（略） 4，数据分析 max(x)求x各列的最大元素 mean(x)求x各列的平均值 median(x)找出x各列的中位元素 min(x)求出x各列的最小元素 S=cumsum(x)求x各列元素累计和 sort(x)使x的各列元素按递增排序。 std(x)求x各列的标准差。 sum(x)求x各列元素之和 prod(x)求x各列元素之积
第6章 MATLAB 数值计算 例6.1 求矩阵A 的每行及每列的最大和最小元素，并求整个矩阵的最大和最小元素。 1356 78256323578255631 01-???? -? ?=???? -??A A=[13,-56,78;25,63,-235;78,25,563;1,0,-1]; max(A,[],2) %求每行最大元素 min(A,[],2) %求每行最小元素 max(A) %求每列最大元素 min(A) %求每列最小元素 max(max(A)) %求整个矩阵的最大元素。也可使用命令：max(A(:)) min(min(A)) %求整个矩阵的最小元素。也可使用命令：min(A(:)) 例6.2 求矩阵A 的每行元素的乘积和全部元素的乘积。 A=[1,2,3,4;5,6,7,8;9,10,11,12]; S=prod(A,2) prod(S) %求A 的全部元素的乘积。也可以使用命令prod(A(:)) 例6.3 求向量X =(1！，2！，3！，…，10！)。 X=cumprod(1:10) 例6.4 对二维矩阵x ，从不同维方向求出其标准方差。 x=[4,5,6;1,4,8] %产生一个二维矩阵x y1=std(x,0,1) y2=std(x,1,1) y3=std(x,0,2) y4=std(x,1,2) 例6.5 生成满足正态分布的10000×5随机矩阵，然后求各列元素的均值和标准方差，再求这5列随机数据的相关系数矩阵。 X=randn(10000,5); M=mean(X) D=std(X) R=corrcoef(X)
例6.6 对下列矩阵做各种排序。 185412613713-?? ??=?? ??-?? A A=[1,-8,5;4,12,6;13,7,-13]; sort(A) %对A 的每列按升序排序 -sort(-A,2) %对A 的每行按降序排序 [X,I]=sort(A) %对A 按列排序，并将每个元素所在行号送矩阵I 例6.7 给出概率积分 2 (d x x f x x -? e 的数据表如表6.1所示，用不同的插值方法计算f (0.472)。 x=0.46:0.01:0.49; %给出x ，f(x) f=[0.4846555,0.4937542,0.5027498,0.5116683]; format long interp1(x,f,0.472) %用默认方法，即线性插值方法计算f(x) interp1(x,f,0.472,'nearest') %用最近点插值方法计算f(x) interp1(x,f,0.472,'spline') %用3次样条插值方法计算f(x) interp1(x,f,0.472,'cubic') %用3次多项式插值方法计算f(x) format short 例6.8 某检测参数f 随时间t 的采样结果如表6.2，用数据插值法计算t =2，7，12，17，22，17，32，37，42，47，52，57时的f 值。 T=0:5:65; X=2:5:57;
列主元法 function lianzhuyuan(A,b) n=input('请输入n:') %选择阶数A=zeros(n,n); %系数矩阵A b=zeros(n,1); %矩阵b X=zeros(n,1); %解X for i=1:n for j=1:n A(i,j)=(1/(i+j-1)); %生成hilbert矩阵A end b(i,1)=sum(A(i,:)); %生成矩阵b end for i=1:n-1 j=i; top=max(abs(A(i:n,j))); %列主元 k=j; while abs(A(k,j))~=top %列主元所在行 k=k+1; end for z=1:n %交换主元所在行a1=A(i,z); A(i,z)=A(k,z); A(k,z)=a1; end a2=b(i,1); b(i,1)=b(k,1); b(k,1)=a2; for s=i+1:n %消去算法开始m=A(s,j)/A(i,j); %化简为上三角矩阵 A(s,j)=0; for p=i+1:n A(s,p)=A(s,p)-m*A(i,p); end b(s,1)=b(s,1)-m*b(i,1); end end X(n,1)=b(n,1)/A(n,n); %回代开始 for i=n-1:-1:1 s=0; %初始化s for j=i+1:n s=s+A(i,j)*X(j,1);
end X(i,1)=(b(i,1)-s)/A(i,i); end X 欧拉法 clc clear % 欧拉法 p=10; %贝塔的取值 T=10; %t取值的上限 y1=1; %y1的初值 r1=1; %y2的初值 %输入步长h的值 h=input('欧拉法please input number(h=1 0.5 0.25 0.125 0.0625):h=') ; if h>1 or h<0 break end S1=0:T/h; S2=0:T/h; S3=0:T/h; S4=0:T/h; i=1; % 迭代过程 for t=0:h:T Y=(exp(-t)); R=(1/(p-1))*exp(-t)+((p-2)/(p-1))*exp(-p*t); y=y1+h*(-y1); y1=y; r=r1+h*(y1-p*r1); r1=r; S1(i)=Y; S2(i)=R; S3(i)=y; S4(i)=r; i=i+1; end t=[0:h:T]; % 红线为解析解，'x'为数值解 plot(t,S1,'r',t,S3,'x')
Matlab关于数值计算的实现 摘要：数值计算（numerical computation computation），主要研究更好的利用计算机更好的进行数值计算，解决各种数学问题。数值分析包括离散傅里叶变换，考虑截断误差，计算误差，函数的敛散性与稳定性等。在数学方面，数值计算的主要研究数值微分与积分，数据的处理与多项式计算，最优化问题，线性方程与非线性方程的求解，常微分方程的数值求解等。同时，数值计算在物理，化学，经济等方面也有研究，本文暂且不表。M atlab软件历经二十多年来的发展，已成为风靡世界的数学三大软件（matlb，Mathematica l，Maple）之一，在数学类科技应用软件中在数值计算方面首屈一指。Matlab以矩阵为数据操作的基本单位，使得矩阵运算十分便捷快速，同时Matlab还提供了海量的计算函数，而且使用可靠地算法进行计算，能使用户在繁复的数学运算中解脱，Matlab还具有方便且完善的图形处理功能，方便绘制二维和三维图形并修饰。
目录 1.数值分析(离散傅里叶变换，考虑截断误差，计算误差，函数 的敛散性与稳定性) 2.数值计算(数值微分与积分，数据的处理与多项式计算， 最优化问题，线性方程与非线性方程的求解，常微分方程的数值求解) 3.图形处理功能(方便绘制二维和三维图形并修饰) 4.总结
1.数据统计与分析 Matlab 可以进行求矩阵的最大最小元素，平均值与中值，关于矩阵元素的求和与求积，累加和与累乘积，标准方程，相关系数，元素排序。现在以求标准方差举例说明Matlab 的实现。 在Matlab 中，实现标准方差计算的函数为std 。对于向量(Y ），std （Y ）实现返回一个标准方差，而对于矩阵（A ），std （A ）返回一个行向量，该行向量的每个元素对应着矩阵A 各行或各列的标准方差。一般调用std 函数的格式为std （A ，flag ，dim ） Dim 取1或者2分别对应求各列或各行的标准方差，flag 取1时，按照标准方差的计算公式 ∑-=-=N i x x S i N 1 2 1)(11来计算。若flag 取2，则用公式 ∑-==N i x x S i N 1 2 2) (1 进行计算。默认的flag 取值为0，dim 取值为1。课本page143 2. 离散傅里叶变换 离散傅里叶变换广泛应用于信号的分析，光谱和声谱分析、全息技术等各个领域。但直接计算dft 的运算量与变化的长度N 的平方成正比，当N 较大时，计算量太大。随着计算机技术的迅速发展，在计算机上进行离散傅里叶变换计算成为可能。特别是快速傅里叶变换算法的出现，为傅里叶变换的应创造了条件。 （1）：傅里叶变换算法的简述。 傅立叶变换是一种分析信号的方法，它可分析信号的成分，也可用这些成分合成信号。许多波形可作为信号的成分，比如正弦波、方波、锯齿波等，傅立叶变换用正弦波作为信号的成分. f(t)是t 的周期函数，如果t 满足狄里赫莱条件：在一个以2T 为周期内f(X)连续或只有有限个第一类间断点，附f （x ）单调或可划分成有限个单调区间，则F （x ）以2T 为周期的傅里叶级数收敛，和函数S （x ）也是以2T 为周期的周期函数，且在这些间断点上，函数是有限值；在一个周期内具有有限个极值点；绝对可积。则有下图①式成立。称为积分运算f(t)的傅立叶变换， ②式的积分运算叫做F(ω)的傅立叶逆变换。F(ω)叫做f(t)的像函数，f(t)叫做 F(ω)的像原函数。F(ω)是f(t)的像。f(t)是F(ω)原像。 ①傅立叶变换 ②傅立叶逆变换
课程名称：Matlab语言 开设时间：2016—2017学年第 2 学期 专业班级：学生学号：学生姓名： 实验名称：实验二、MATLAB的数值运算和程序实验成绩： 指导教师：批改时间： 一、实验目的和要求 1）掌握基本的矩阵运算及常用的函数。 2）掌握MATLAB函数的编写及调试方法。 3）掌握MATLAB常用的数值运算函数。 二、实验仪器和设备 计算机一台 三、实验过程 1、一维数组在命令窗口执行下面指令，观察输出结果，体味数组创建和寻访方法，%号后面的为注释，不用输入。 rand('state',0) % 把均匀分布伪随机发生器置为0 状态 x=rand(1,5) % 产生（1*5）的均布随机数组 x(3) % 寻访数组x 的第三个元素。 x([1 2 5]) % 寻访数组x 的第一、二、五个元素组成的子数组。 x(1:3) % 寻访前三个元素组成的子数组 x(3:end) % 寻访除前2 个元素外的全部其他元素。end 是最后一个元素的下标。 x(3:-1:1) % 由前三个元素倒排构成的子数组 x(find(x>0.5)) % 由大于0.5 的元素构成的子数组 x([1 2 3 4 4 3 2 1]) % 对元素可以重复寻访，使所得数组长度允许大于原数组。 x(3) = 0 % 把上例中的第三个元素重新赋值为0。 x([1 4])=[1 1] % 把当前x 数组的第一、四个元素都赋值为1。 x[3]=[] % 空数组的赋值操作
2、在命令窗口执行下面指令，观察输出结果 a=2.7358; b=33/79; % 这两条指令分别给变量 a ， b 赋值。 C=[1,2*a+i*b,b*sqrt(a);sin(pi/4),a+5*b,3.5+i] % 这指令用于创建二维组C M_r=[1,2,3;4,5,6],M_i=[11,12,13;14,15,16] % 创建复数数组的另一种方法 CN=M_r+i*M_i % 由实部、虚部数组构成复数数组 3. 记录下面题目的程序和运行后的结果。 1??????=654321a ??????-=531142b ???? ? ?????-=201c ??????????=063258741d 下列运算是否合法，为什么？如合法，结果是多少？
本文档包含上一个文档中的五个数值分析实验题C语言程序及Matlab程序实验一 C程序 #include "stdio.h" #include "math.h" void main() { inti=0; float a=0.1,b=1.9,t=0.0,e=1.9; if((pow(a,7)-28*pow(a,4)+14)*(pow(b,7)-28*pow(b,4)+14)<0) if((7*pow(x,6)-112*pow(x,3))) printf("x=%f,i=%d,e=%f\n",x,i,e); for(i=1;i<7&&e>0.00001;i++) { t=x; x=x-(pow(x,7)-28*pow(x,4)+14)/(7*pow(x,6)-112*pow(x,3)); e=fabs(t-x); printf("x=%f,i=%d,e=%f\n",x,i,e);
} } Matable 程序 i=0; x=1.9;t=0.0;e=1.9; disp(['i=',num2str(i),' ','x=',num2str(x),' ','e=',num2str(e)]); for i=1:7 t=x; x=x-(x^7-28*x^4+14)/(7*x^6-112*x^3); e=abs(t-x); disp(['i=',num2str(i),' ','x=',num2str(x),' ','e=',num2str(e)]); if e<0.00001 break; end end 实验二 C程序 #include"stdio.h" #include"math.h" //已知量 double x[10]={1,2,3,4,5,6,7,8,9,10};
function [x0,k]=bisect1(fun1,a,b,ep) if nargin<4 ep=1e-5; end fa=feval(fun1,a); fb=feval(fun1,b); if fa*fb>0 x0=[fa,fb]; k=0; return; end k=1; while abs(b-a)/2>ep x=(a+b)/2; fx=feval(fun1,x); if fx*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) if nargin<4 N=500; end if nargin<3 ep=1e-5; end x=x0; x0=x+2*ep;
while abs(x-x0)>ep & k 数值分析重要算法的matlab程序 插值多项式:拉格朗日插值 function yh=lagrange(x,y,xh) n=length(x); m=length(xh); yh=zeros(1,m); c1=ones(n-1,1); c2=ones(1,m); for i=1:n xp=x([1:i-1 i+1:n]); yh=yh+y(i)*prod((c1*xh-xp'*c2)./(x(i)-xp'*c2)); end 插值多项式:牛顿插值（以数值实验3.2）为例 function y=ex32 n = 21; x = linspace(-5,5,n)'; h = (5-(-5))/(n-1); y = 1./(1+x.^2); % form the differences table for j = 2:n, y(1:n+1-j,j) = diff(y(1:n+2-j,j-1))./(x(j:n)-x(1:n+1-j)); end % newton coeff y = y(1,:); pz = [ ]; v = linspace(-5,5,80); for t = v, z = y(n); for j = n-1:-1:1, z = z * (t-x(j))+y(j); end pz=[pz z]; end plot(v,pz,'r+-',v,1./(1+v.^2),'g--'); 数值积分：梯形求积公式求积分 function I=ftrapz(fun,a,b,n) h=(b-a)/n; x=linspace(a,b,n+1); y=feval(fun,x); I=h*(0.5*y(1)+sum(y(2:n))+0.5*y(n+1)); 数值积分：抛物型求积公式求积分 function I=fsimpsion(fun,a,b,n) h=(b-a)/n; x=linspace(a,b,2*n+1); y=feval(fun,x); I=(h/6)*(y(1)+2*sum(y(3:2:2*n-1))+4*sum(y(2:2:2*n))+y(2*n+1)); 追赶法解三对角方程组 Gauss_Seidel迭代法的Matlab程序 function [x]=Gauss_Seidel_iterative(A,b) % 用Gauss_Seidel迭代求解线性方程组,矩阵A是方阵 x0=zeros(1,length(b)); % 赋初值 tol=10^(-2); % 给定误差界 N=1000; % 给定最大迭代次数 [n,n]=size(A); % 确定矩阵A的阶 k=1; % 迭代过程 while k<=N x(1)=(b(1)-A(1,2:n)*x0(2:n)')/A(1,1); for i=2:n x(i)=(b(i)-A(i,1:i-1)*x(1:i-1)'-A(i,i+1:n)*x0(i+1:n)')... /A(i,i); end if max(abs(x-x0))<=tol fid = fopen('G_S_iter_result.txt', 'wt'); fprintf(fid,'\n********用Gauss_Seidel迭代求解线性方程组的输出结果********\n\n'); fprintf(fid,'迭代次数: %d次\n\n',k); fprintf(fid,'x的值\n\n'); fprintf(fid, '%12.8f \n', x); break; end k=k+1; x0=x; end if k==N+1 fid = fopen('G_S_iter_result.txt', 'wt'); fprintf(fid,'\n********用Gauss_Seidel迭代求解线性方程组的输出结果********\n\n'); fprintf(fid,'迭代次数: %d次\n\n',k); 第四章数值计算 4.1引言 本章将花较大的篇幅讨论若干常见数值计算问题：线性分析、一元和多元函数分析、微积分、数据分析、以及常微分方程（初值和边值问题）求解等。但与一般数值计算教科书不同，本章的讨论重点是：如何利用现有的世界顶级数值计算资源MATLAB。至于数学描述，本章将遵循“最低限度自封闭”的原则处理，以最简明的方式阐述理论数学、数值数学和MATLAB计算指令之间的在联系及区别。 对于那些熟悉其他高级语言（如FORTRAN，Pascal，C++）的读者来说，通过本章，MATLAB 卓越的数组处理能力、浩瀚而灵活的M函数指令、丰富而友善的图形显示指令将使他们体验到解题视野的豁然开朗，感受到摆脱烦琐编程后的眉眼舒展。 对于那些经过大学基本数学教程的读者来说，通过本章，MATLAB精良完善的计算指令，自然易读的程序将使他们感悟“教程”数学的基础地位和局限性，看到从“理想化”简单算例通向科学研究和工程设计实际问题的一条途径。 对于那些熟悉MATLAB基本指令的读者来说，通过本章，围绕基本数值问题展开的容将使他们体会到各别指令的运用场合和在关系，获得综合运用不同指令解决具体问题的思路和借鉴。 由于MATLAB的基本运算单元是数组，所以本章容将从矩阵分析、线性代数的数值计算开始。然后再介绍函数零点、极值的求取，数值微积分，数理统计和分析，拟合和插值，Fourier分析，和一般常微分方程初值、边值问题。本章的最后讨论稀疏矩阵的处理，因为这只有在大型问题中，才须特别处理。 从总体上讲，本章各节之间没有依从关系，即读者没有必要从头到尾系统阅读本章容。读者完全可以根据需要阅读有关节次。除特别说明外，每节中的例题指令是独立完整的，因此读者可以很容易地在自己机器上实践。 MATLAB从5.3版升级到6.x版后，本章容的变化如下： ●MATLAB从6.0版起，其矩阵和特征值计算指令不再以LINPACK和EISPACK库为基础， 而建筑在计算速度更快、运行更可靠的LAPACK和ARPACK程序库的新基础上。因此，虽然各种矩阵计算指令没有变化，但计算结果却可能有某些不同。这尤其突出地表现在涉及矩阵分解、特征向量、奇异向量等的计算结果上。对此，用户不必诧异，因为构成空间的基向量时不唯一的，且新版的更可信。本书新版全部算例结果是在6.x版上给出的。 ●在5.3版本中，泛函指令对被处理函数的调用是借助函数名字符串进行的。这种调用 方式在6.x版中已被宣布为“过渡期允许使用但即将被淘汰的调用方式”；而新的调用方式是借助“函数句柄”进行的。因此，关于述泛函指令，本章新版着重讲述如何使用“函数句柄”，同时兼顾“函数名字符串”调用法。 ●MATLAB从6.0版起，提供了一组专门求微分方程“边值问题”数值解的指令。适应这 种变化，本章新增第4.14.5节，用2个算例阐述求解细节。 ● 5.3版中的积分指令quad8已经废止；6.x版启用新积分指令quad l；6.5版新增三重 积分指令triplequad。本章新版对此作了相应的改变。 4.2LU分解和恰定方程组的解 1，秦九韶算法，求出P(x=3)=2+4x+5x^2+2x^3的值 clear all; x=3; n=3; a(1)=2;a(2)=4;a(3)=5;a(4)=2 v(1)=a(n+1); for k=2:(n+1); v(k)=x*v(k-1)+a(n-k+2); end p=v(n+1) p =, 113 2，一次线型插值程序：利用100.121.求115的开方。 clear all; x1=100;x2=121;y1=10;y2=11; x=115; l1=(x-x2)/(x1-x2);l2=(x-x1)/(x2-x1); p1=l1*y1+l2*y2 p1 = 10.7143 3，分段插值程序，已知为S1（x)为（0，0），（1，1），（2，5）（3, 8）上的分段一次插值，求S1（1.5）. clear all x=[0 1 2 3]; y=[0 1 5 8]; n=length(x);a=1.5; for i=2:n if(x(i-1)<=a end H1=y(i-1)+(y(i)-y(i-1))/(x(i)-x(i-1))*(a-x(i-1)) H1 = 3.5000 4）曲线拟合：用一个5次多项式在区间[0，2π]内逼近函数sin(x)。clear all X=linspace(0,2*pi,50);Y=sin(X); [P,S]=polyfit(X,Y,5) plot(X,Y,'k*',X,polyval(P,X),'k-') P = -0.0056 0.0874 -0.3946 0.2685 0.8797 0.0102 S = R: [6x6 double] df: 44 normr: 0.0337 5)求有理分式的导数 clear all P=[3,5,0,-8,1,-5]; Q=[10,5,0,0,6,0,0,7,-1,0,-100]; [p,q]=polyder(P,Q) 6）将以下数据按从小到大排序：4.3 5.7 5.2 1.8 9.4 a=[4.3 5.7 5.2 1.8 9.4];b(1:100)=0;n=1; b(a*10)=1; for k=1:100 a(n)=k/10; if b(k)>0 a(n)=k/10; n=n+1; end end a a = 2012 年数学建模培训材料——Matlab 软件的使用 第一讲Matlab基本数值计算 一、矩阵 在Matlab 中，一个矩阵可以使数学意义上的矩阵，也可以是标量或者向量。对于一个标量（一个数）可以将之作为 1 1的矩阵，而向量（一行或一列）则可以认为是 1 n 或者 n 1的矩阵。另外，一个 0 0 矩阵在 Matlab 中被认为是空矩阵，用“[] ”表示。 1、矩阵的创建 矩阵的创建可以有以下几种形式 ⑴直接输入 >> A=[1 2 3;4 3 7;2 4 1] 注意：每行间的元素用逗号或空格分开，行与行之间用分号或回车分开，矩阵标示是一对中括号[ ] 。 也可以采用数组编辑器（Array Editor ）像在 Excel 电子表格中据那样输入数据。 ⑵通过语句和函数产生 常用的特殊矩阵： zeros:全零矩阵 ,ones：全 1 矩阵， eye：单位矩阵， rand:随机矩阵， diag:对角阵等。 例： >> A=ones(3,4) >>E=eye(3) >>D=diag([3 5 2]) ⑶对矩阵进行裁剪或拼接 ⑷从外部文件装入数据 外部数据文件可以是以保存的Matlab 工作空间，也可以是文本（.txt）文件，或者是电子表格创建的文件（ .xls）. 例：已知一个文本格式的数据文件 E:\Mathmodel\data1.txt >> load e:\Mathmodel\data1.txt 得到一个变量名与文件名相同的矩阵（data1）。注意：文件的扩展名不能省略。 例：已知一个 Excel 文件的路径为E:\Mathmodel\data2.xls a. 缺省操作 : >> NUMBER=xlsread('E:\Mathmodel\data2.xls') >>[NUMBER,TXT]=xlsread('E:\Mathmodel\data2.xls') 默认操作是从第一个工作表（sheet1）中提取数据。 b.从指定的工作表（而不是第一个）中提取数据： >> NUMBER=xlsread('E:\Mathmodel\data2.xls','S2') 或者 >>NUMBER=xlsread('E:\Mathmodel\data2.xls',2) c.从指定的工作表中读取指定区域的数据： >>NUMBER=xlsread('E:\Mathmodel\data2.xls',2,'g3:i8') 2、Matlab 的矩阵运算 ⑴基本运算 矩阵的加（ +）、减（ -）、乘（ * ）、乘方（ ^）运算法则与代数中的定义完全一致。例如： >>A=[1 2;3 4];B=[3 1;4 8];数值分析重要算法的matlab程序
matlab数值计算(diedai)
matlab入门经典教程--第四章 数值计算
数值分析matlab程序实例
(完整版)第一讲Matlab基本数值计算.doc
相关主题
文本预览