MATLAB在数值分析中的应用

MATLAB在数值分析中的应用
MATLAB在数值分析中的应用

MATLAB 在数值分析中的应用

插值与拟合是来源于实际、又广泛应用于实际的两种重要方法。随着计算机的不断发展及计算水平的不断提高,它们已在国民生产和科学研究等方面扮演着越来越重要的角色。下面对插值中分段线性插值、拟合中的最为重要的最小二乘法拟合加以介绍。

1 分段线性插值

所谓分段线性插值就是通过插值点用折线段连接起来逼近原曲线,这也是计算机绘制图形的基本原理。实现分段线性插值不需编制函数程序,MA TLAB 自身提供了内部函数interp1其主要用法如下:

interp1(x,y,xi) 一维插值

◆ yi=interp1(x,y,xi)

对一组点(x,y) 进行插值,计算插值点xi 的函数值。x 为节点向量值,y 为对应的节点函数值。如果y 为矩阵,则插值对y 的每一列进行,若y 的维数超出x 或 xi 的维数,则返回NaN 。

◆ yi=interp1(y,xi)

此格式默认x=1:n ,n 为向量y 的元素个数值,或等于矩阵y 的size(y,1)。

◆ yi=interp1(x,y,xi,’method’)

method 用来指定插值的算法。默认为线性算法。其值常用的可以是如下的字符串。 ● nearest 线性最近项插值。

● linear 线性插值。

● spline 三次样条插值。

● cubic 三次插值。

所有的插值方法要求x 是单调的。x 也可能并非连续等距的。

正弦曲线的插值示例:

>> x=0:0.1:10;

>> y=sin(x);

>> xi=0:0.25:10;

>> yi=interp1(x,y,xi);

>> plot(x,y,’0’,xi,yi)

则可以得到相应的插值曲线(读者可自己上机实验)。

Matlab 也能够完成二维插值的运算,相应的函数为interp2,使用方法与interpl 基本相同,只是输入和输出的参数为矩阵,对应于二维平面上的数据点,详细的用法见Matlab 联机帮助。

2 最小二乘法拟合

在科学实验的统计方法研究中,往往要从一组实验数据),(i i y x 中寻找出自变量x 和因

变量y 之间的函数关系y=f(x) 。由于观测数据往往不够准确,因此并不要求y=f(x)经过所有的点 ),(i i y x ,而只要求在给定点i x 上误差i i i y x f -=)(δ按照某种标准达到最小,通常采用欧氏范数2δ作为误差量度的标准。这就是所谓的最小二乘法。在MA TLAB 中实现最小二乘法拟合通常采用polyfit 函数进行。

函数polyfit 是指用一个多项式函数来对已知数据进行拟合,我们以下列数据为例介绍

这个函数的用法:

>> x=0:0.1:1;

>> y=[ -0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2 ]

为了使用polyfit,首先必须指定我们希望以多少阶多项式对以上数据进行拟合,如果我们指定一阶多项式,结果为线性近似,通常称为线性回归。我们选择二阶多项式进行拟合。

>> P= polyfit (x, y, 2)

P=

-9.8108 20.1293 -0.0317

函数返回的是一个多项式系数的行向量,写成多项式形式为:

.92-

+

x

-x

8108

0317

20

.0

1293

.

为了比较拟合结果,我们绘制两者的图形:

>> xi=linspace (0, 1, 100); %绘图的X-轴数据。

>> Z=polyval (p, xi); %得到多项式在数据点处的值。

当然,我们也可以选择更高幂次的多项式进行拟合,如10阶:

>> p=polyfit (x, y, 10);

>> xi=linspace (0, 1,100);

>> z=ployval (p, xi);

读者可以上机绘图进行比较,曲线在数据点附近更加接近数据点的测量值了,但从整体上来说,曲线波动比较大,并不一定适合实际使用的需要,所以在进行高阶曲线拟合时,“越高越好”的观点不一定对的。

《MATLAB与数值分析》第一次上机实验报告

电子科技大学电子工程学院标准实验报告(实验)课程名称MATLAB与数值分析 学生姓名:李培睿 学号:2013020904026 指导教师:程建

一、实验名称 《MATLAB与数值分析》第一次上机实验 二、实验目的 1. 熟练掌握矩阵的生成、加、减、乘、除、转置、行列式、逆、范数等运算 操作。(用.m文件和Matlab函数编写一个对给定矩阵进行运算操作的程序) 2. 熟练掌握算术符号操作和基本运算操作,包括矩阵合并、向量合并、符号 转换、展开符号表达式、符号因式分解、符号表达式的化简、代数方程的符号解析解、特征多项式、函数的反函数、函数计算器、微积分、常微分方程的符号解、符号函数的画图等。(用.m文件编写进行符号因式分解和函数求反的程序) 3. 掌握Matlab函数的编写规范。 4、掌握Matlab常用的绘图处理操作,包括:基本平面图、图形注释命令、 三维曲线和面的填充、三维等高线等。(用.m文件编写在一个图形窗口上绘制正弦和余弦函数的图形,并给出充分的图形注释) 5. 熟练操作MATLAB软件平台,能利用M文件完成MATLAB的程序设计。 三、实验内容 1. 编程实现以下数列的图像,用户能输入不同的初始值以及系数。并以x, y为坐标显示图像 x(n+1) = a*x(n)-b*(y(n)-x(n)^2); y(n+1) = b*x(n)+a*(y(n)-x(n)^2) 2. 编程实现奥运5环图,允许用户输入环的直径。 3. 实现对输入任意长度向量元素的冒泡排序的升序排列。不允许使用sort 函数。 四、实验数据及结果分析 题目一: ①在Editor窗口编写函数代码如下:

数值分析MATLAB上机实验

数值分析实习报告 姓名:gestepoA 学号:201******* 班级:***班

序言 随着计算机技术的迅速发展,数值分析在工程技术领域中的应用越来越广泛,并且成为数学与计算机之间的桥梁。要解决工程问题,往往需要处理很多数学模型,不仅要研究各种数学问题的数值解法,同时也要分析所用的数值解法在理论上的合理性,如解法所产生的误差能否满足精度要求:解法是否稳定、是否收敛及熟练的速度等。而且还能减少大量的人工计算。 由于工程实际中所遇到的数学模型求解过程迭代次数很多,计算量很大,所以需要借助如MATLAB,C++,VB,JAVA的辅助软件来解决,得到一个满足误差限的解。本文所计算题目,均采用MATLAB进行编程,MATLAB被称为第四代计算机语言,利用其丰富的函数资源,使编程人员从繁琐的程序代码中解放出来MATLAB最突出的特点就是简洁,它用更直观的、符合人们思维习惯的代码。它具有以下优点: 1友好的工作平台和编程环境。MATLAB界面精致,人机交互性强,操作简单。 2简单易用的程序语言。MATLAB是一个高级的矩阵/阵列语言,包含控制语言、函数、数据结构,具有输入、输出和面向对象编程特点。用户可以在命令窗口中将输入语句与执行命令同步,也可以先编好一个较大的复杂的应用程序(M 文件)后再一起运行。 3强大的科学计算机数据处理能力。包含大量计算算法的集合,拥有600多个工程中要用到的数学运算函数。 4出色的图像处理功能,可以方便地输出二维图像,便于我们绘制函数图像。

目录 1 第一题 (4) 1.1 实验目的 (4) 1.2 实验原理和方法 (4) 1.3 实验结果 (5) 1.3.1 最佳平方逼近法 (5) 1.3.2 拉格朗日插值法 (7) 1.3.3 对比 (8) 2 第二题 (9) 2.1实验目的 (9) 2.2 实验原理和方法 (10) 2.3 实验结果 (10) 2.3.1 第一问 (10) 2.3.2 第二问 (11) 2.3.3 第三问 (11) 3 第三题 (12) 3.1实验目的 (12) 3.2 实验原理和方法 (12) 3.3 实验结果 (12) 4 MATLAB程序 (14)

matlab数值计算(命令与示例)

MATLAB数值计算 MATLAB数值计算 (1) 1创建矩阵 (3) 1.1直接输入 (3) 1.2向量 (3) 1.2.1linspace:线性分布 (3) 1.2.2冒号法 (3) 1.3函数创建 (4) 1.3.1eye:单位矩阵 (4) 1.3.2rand:随机矩阵 (4)

1.3.3zeros:全0矩阵 (4) 1.3.4ones:全1矩阵 (5) 2矩阵运算 (5) 2.1加减 (5) 2.1.1[M×N]±[M×N] (5) 2.2乘 (6) 2.2.1[M×N]*a (6) 2.2.2[M×N]*[N×M] (6) 2.3乘方 (7) 2.3.1[M×M]^a (7) 2.3.2a^[M×M] (7) 2.4特殊运算 (8) 2.4.1求逆inv (8) 2.4.2行列式det (8) 2.4.3特征值eig (8) 2.4.4转置'和.' (9) 2.4.5变形reshape (10) 2.4.6翻转rot90,fliplr,flipud (11) 2.4.7抽取diag,tril,triu (12) 2.5数组运算 (12) 2.5.1乘 (12) [M×N].*[M×N] (12) 2.5.2除 (13) [M×N]./[M×N] (14) [M×N].\[M×N] (14) 2.5.3乘方 (14) [M×N].^[M×N] (15) a.^[M×N] (15) 2.6除法 (15) 2.6.1求解线性方程组 (15) 3多项式 (16) 3.1系数表示法poly (16) 3.2求根roots (16) 3.3乘法conv (16) 3.4除法deconv (17) 3.5求值polyval (17) 3.6微分polyder (18)

数值分析的matlab实现

第2章牛顿插值法实现 参考文献:[1]岑宝俊. 牛顿插值法在凸轮曲线修正设计中的应用[J]. 机械工程师,2009,10:54-55. 求牛顿插值多项式和差商的MA TLAB 主程序: function[A,C,L,wcgs,Cw]=newpoly(X,Y) n=length(X);A=zeros(n,n);A(:,1) =Y'; s=0.0;p=1.0;q=1.0;c1=1.0; for j=2:n for i=j:n A(i,j)=(A(i,j-1)-A(i-1,j-1))/(X(i)-X(i-j+1)); end b=poly(X(j-1));q1=conv(q,b);c1=c1*j;q=q1; end C=A(n,n);b=poly(X(n));q1=conv(q1,b); for k=(n-1):-1:1 C=conv(C,poly(X(k)));d=length(C);C(d)=C(d)+A(k,k); end L(k,:)=poly2sym(C);Q=poly2sym(q1); syms M wcgs=M*Q/c1;Cw=q1/c1; (1)保存名为newpoly.m 的M 文件 (2)输入MA TLAB 程序 >> X=[242,243,249,250]; >> Y=[13.681,13.526,13.098,13.095]; >> [A,C,L,wcgs,Cw]=newpoly(X,Y) 输出3阶牛顿插值多项式L 及其系数向量C 差商的矩阵A ,插值余项wcgs 及其 ) ()()1(ξ+n n f x R 的系数向量Cw 。 A = 13.6810 0 0 0 13.5260 -0.1550 0 0 13.0980 -0.0713 0.0120 0 13.0950 -0.0030 0.0098 -0.0003 C = 1.0e+003 *

Matlab作业3(数值分析)答案

Matlab作业3(数值分析) 机电工程学院(院、系)专业班组 学号姓名实验日期教师评定 1.计算多项式乘法(x2+2x+2)(x2+5x+4)。 答: 2. (1)将(x-6)(x-3)(x-8)展开为系数多项式的形式。(2)求解在x=8时多项 式(x-1)(x-2) (x-3)(x-4)的值。 答:(1) (2)

3. y=sin(x),x从0到2π,?x=0.02π,求y的最大值、最小值、均值和标准差。 4.设x=[0.00.30.8 1.1 1.6 2.3]',y=[0.500.82 1.14 1.25 1.35 1.40]',试求二次多项式拟合系数,并据此计算x1=[0.9 1.2]时对应的y1。解:x=[0.0 0.3 0.8 1.1 1.6 2.3]'; %输入变量数据x y=[0.50 0.82 1.14 1.25 1.35 1.40]'; %输入变量数据y p=polyfit(x,y,2) %对x,y用二次多项式拟合,得到系数p x1=[0.9 1.2]; %输入点x1 y1=polyval(p,x1) %估计x1处对应的y1 p = -0.2387 0.9191 0.5318 y1 = a) 1.2909

5.实验数据处理:已知某压力传感器的测试数据如下表 p为压力值,u为电压值,试用多项式 d cp bp ap p u+ + + =2 3 ) ( 来拟 合其特性函数,求出a,b,c,d,并把拟合曲线和各个测试数据点画在同一幅图上。解: >> p=[0.0,1.1,2.1,2.8,4.2,5.0,6.1,6.9,8.1,9.0,9.9]; u=[10,11,13,14,17,18,22,24,29,34,39]; x=polyfit(p,u,3) %得多项式系数 t=linspace(0,10,100); y=polyval(x,t); %求多项式得值 plot(p,u,'*',t,y,'r') %画拟和曲线 x = 0.0195 -0.0412 1.4469 9.8267

matlab与数值分析作业

数值分析作业(1) 1:思考题(判断是否正确并阐述理由) (a)一个问题的病态性如何,与求解它的算法有关系。 (b)无论问题是否病态,好的算法都会得到它好的近似解。 (c)计算中使用更高的精度,可以改善问题的病态性。 (d)用一个稳定的算法计算一个良态问题,一定会得到他好的近似解。 (e)浮点数在整个数轴上是均匀分布。 (f)浮点数的加法满足结合律。 (g)浮点数的加法满足交换律。 (h)浮点数构成有效集合。 (i)用一个收敛的算法计算一个良态问题,一定得到它好的近似解。√2: 解释下面Matlab程序的输出结果 t=0.1; n=1:10; e=n/10-n*t 3:对二次代数方程的求解问题 20 ++= ax bx c 有两种等价的一元二次方程求解公式

2224b x a c x b ac -±==- 对 a=1,b=-100000000,c=1,应采用哪种算法? 4:函数sin x 的幂级数展开为: 357 sin 3!5!7! x x x x x =-+-+ 利用该公式的Matlab 程序为 function y=powersin(x) % powersin. Power series for sin(x) % powersin(x) tries to compute sin(x)from a power series s=0; t=x; n=1; while s+t~=s; s=s+t; t=-x^2/((n+1)*(n+2))*t n=n+2; end

(a ) 解释上述程序的终止准则; (b ) 对于x=/2π、x=11/2π、x =21/2π,计算的精度是多少?分别需 要计算多少项? 5:指数函数的幂级数展开 2312!3!x x x e x =+++ + 根据该展开式,编写Matlab 程序计算指数函数的值,并分析计算结果(重点分析0x <的计算结果)。

数值分析的MATLAB程序

列主元法 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与数值分析课程总结

MATLAB与数值分析课程总结 姓名:董建伟 学号:2015020904027 一:MATLAB部分 1.处理矩阵-容易 矩阵的创建 (1)直接创建注意 a中括号里可以用空格或者逗号将矩阵元素分开 b矩阵元素可以是任何MATLAB表达式,如实数复数等 c可以调用赋值过的任何变量,变量名不要重复,否则会被覆盖 (2)用MATLAB函数创建矩阵如:a空阵[] b rand/randn——随机矩阵 c eye——单位矩阵 d zeros ——0矩阵 e ones——1矩阵 f magic——产生n阶幻方矩阵等 向量的生成 (1)用冒号生成向量 (2)使用linspace和logspace分别生成线性等分向量和对 数等分向量 矩阵的标识和引用 (1)向量标识 (2)“0 1”逻辑向量或矩阵标识 (3)全下标,单下标,逻辑矩阵方式引用 字符串数组 (1)字符串按行向量进行储存 (2)所有字符串用单引号括起来 (3)直接进行创建 矩阵运算 (1)注意与数组点乘,除与直接乘除的区别,数组为乘方对应元素的幂

(2)左右除时斜杠底部靠近谁谁是分母 (3)其他运算如,inv矩阵求逆,det行列式的值, eig特征值,diag 对角矩阵 2.绘图-轻松 plot-绘制二维曲线 (1)plot(x)绘制以x为纵坐标的二维曲线 plot(x,y) 绘制以x为横坐标,y为纵坐标的二维曲线 x,y为向量或矩阵 (2)plot(x1,y1,x2,y2,。。。。。。)绘制多条曲线,不同字母代替不同颜色:b蓝色,y黄色,r红色,g绿色 (3)hold on后面的pl ot图像叠加在一起 hold off解除hold on命令,plot将先冲去窗口已有图形(4)在hold后面加上figure,可以绘制多幅图形 (5)subplot在同一窗口画多个子图 三维图形的绘制 (1)plot3(x,y,z,’s’) s是指定线型,色彩,数据点形的字 符串 (2)[X,Y]=meshgrid(x,y)生成平面网格点 (3)mesh(x,y,z,c)生成三维网格点,c为颜色矩阵 (4)三维表面处理mesh命令对网格着色,surf对网格片着色 (5)contour绘制二维等高线 (6)axis([x1,xu,y1,yu])定义x,y的显示范围 3.编程-简洁 (1)变量命名时可以由字母,数字,下划线,但是不得包含空格和标点 (2)最常用的数据类型只有双精度型和字符型,其他数据类型只在特殊条件下使用 (3)为得到高效代码,尽量提高代码的向量化程度,避免使用循环结构

数值分析(Hilbert矩阵)病态线性方程组的求解Matlab程序

(Hilbert 矩阵)病态线性方程组的求解 理论分析表明,数值求解病态线性方程组很困难。考虑求解如下的线性方程组的求解Hx = b ,期中H 是Hilbert 矩阵,()ij n n H h ?=,11 ij h i j = +-,i ,j = 1,2,…,n 1. 估计矩阵的2条件数和阶数的关系 2. 对不同的n ,取(1,1,,1)n x =∈K ?,分别用Gauss 消去,Jacobi 迭代,Gauss-seidel 迭 代,SOR 迭代和共轭梯度法求解,比较结果。 3. 结合计算结果,试讨论病态线性方程组的求解。 第1小题: condition.m %第1小题程序 t1=20;%阶数n=20 x1=1:t1; y1=1:t1; for i=1:t1 H=hilb(i); y1(i)=log(cond(H)); end plot(x1,y1); xlabel('阶数n'); ylabel('2-条件数的对数(log(cond(H))'); title('2-条件数的对数(log(cond(H))与阶数n 的关系图'); t2=200;%阶数n=200 x2=1:t2; y2=1:t2; for i=1:t2 H=hilb(i); y2(i)=log(cond(H)); end plot(x2,y2); xlabel('阶数n'); ylabel('2-条件数的对数(log(cond(H))'); title('2-条件数的对数(log(cond(H))与阶数n 的关系图'); 画出Hilbert 矩阵2-条件数的对数和阶数的关系

n=200时 n=20时 从图中可以看出, 1)在n小于等于13之前,图像近似直线 log(cond(H))~1.519n-1.833 2)在n大于13之后,图像趋于平缓,并在一定范围内上下波动,同时随着n的增加稍有上升的趋势 第2小题: solve.m%m第2小题主程序 N=4000;

matlab实现数值分析插值及积分

Matlab实现数值分析插值及积分 摘要: 数值分析(numerical analysis)是研究分析用计算机求解数学计算问题的数值计算方法及其理论的学科,是数学的一个分支,它以数字计算机求解数学问题的理论和方法为研究对象。在实际生产实践中,常常将实际问题转化为数学模型来解决,这个过程就是数学建模。学习数值分析这门课程可以让我们学到很多的数学建模方法。 分别运用matlab数学软件编程来解决插值问题和数值积分问题。题目中的要计算差值和积分,对于问题一,可以分别利用朗格朗日插值公式,牛顿插值公式,埃特金逐次线性插值公式来进行编程求解,具体matlab代码见正文。编程求解出来的结果为:=+。 其中Aitken插值计算的结果图如下: 对于问题二,可以分别利用复化梯形公式,复化的辛卜生公式,复化的柯特斯公式编写程序来进行求解,具体matlab代码见正文。编程求解出来的结果为: 0.6932 其中复化梯形公式计算的结果图如下:

问题重述 问题一:已知列表函数 表格 1 分别用拉格朗日,牛顿,埃特金插值方法计算。 问题二:用复化的梯形公式,复化的辛卜生公式,复化的柯特斯公式计算积分,使精度小于5。 问题解决 问题一:插值方法 对于问题一,用三种差值方法:拉格朗日,牛顿,埃特金差值方法来解决。 一、拉格朗日插值法: 拉格朗日插值多项式如下: 首先构造1+n 个插值节点n x x x ,,,10 上的n 插值基函数,对任一点i x 所对应的插值基函数 )(x l i ,由于在所有),,1,1,,1,0(n i i j x j +-=取零值,因此)(x l i 有因子 )())(()(110n i i x x x x x x x x ----+- 。又因)(x l i 是一个次数不超过n 的多项式,所以只 可能相差一个常数因子,固)(x l i 可表示成: )())(()()(110n i i i x x x x x x x x A x l ----=+- 利用1)(=i i x l 得:

同济大学数值分析matlab编程题汇编

MATLAB 编程题库 1.下面的数据表近似地满足函数2 1cx b ax y ++=,请适当变换成为线性最小二乘问题,编程求最好的系数c b a ,,,并在同一个图上画出所有数据和函数图像. 625 .0718.0801.0823.0802.0687.0606.0356.0995 .0628.0544.0008.0213.0362.0586.0931.0i i y x ---- 解: x=[-0.931 -0.586 -0.362 -0.213 0.008 0.544 0.628 0.995]'; y=[0.356 0.606 0.687 0.802 0.823 0.801 0.718 0.625]'; A=[x ones(8,1) -x.^2.*y]; z=A\y; a=z(1); b=z(2); c=z(3); xh=-1:0.1:1; yh=(a.*xh+b)./(1+c.*xh.^2); plot(x,y,'r+',xh,yh,'b*')

2.若在Matlab工作目录下已经有如下两个函数文件,写一个割线法程序,求出这两个函数 10 的近似根,并写出调用方式: 精度为10 解: >> edit gexianfa.m function [x iter]=gexianfa(f,x0,x1,tol) iter=0; while(norm(x1-x0)>tol) iter=iter+1; x=x1-feval(f,x1).*(x1-x0)./(feval(f,x1)-feval(f,x0)); x0=x1;x1=x; end >> edit f.m function v=f(x) v=x.*log(x)-1; >> edit g.m function z=g(y) z=y.^5+y-1; >> [x1 iter1]=gexianfa('f',1,3,1e-10) x1 = 1.7632 iter1 = 6 >> [x2 iter2]=gexianfa('g',0,1,1e-10) x2 = 0.7549 iter2 = 8

第2讲 matlab的数值分析

第二讲MATLAB的数值分析 2-1矩阵运算与数组运算 矩阵运算和数组运算是MATLAB数值运算的两大类型,矩阵运算是按矩阵的运算规则进行的,而数组运算则是按数组元素逐一进行的。因此,在进行某些运算(如乘、除)时,矩阵运算和数组运算有着较大的差别。在MATLAB中,可以对矩阵进行数组运算,这时是把矩阵视为数组,运算按数组的运算规则。也可以对数组进行矩阵运算,这时是把数组视为矩阵,运算按矩阵的运算规则进行。 1、矩阵加减与数组加减 矩阵加减与数组加减运算效果一致,运算符也相同,可分为两种情况: (1)若参与运算的两矩阵(数组)的维数相同,则加减运算的结果是将两矩阵的对应元素进行加减,如 A=[1 1 1;2 2 2;3 3 3]; B=A; A+B ans= 2 2 2 4 4 4 6 6 6 (2)若参与运算的两矩阵之一为标量(1*1的矩阵),则加减运算的结果是将矩阵(数组)的每一元素与该标量逐一相加减,如 A=[1 1 1;2 2 2;3 3 3]; A+2 ans= 3 3 3 4 4 4 5 5 5 2、矩阵乘与数组乘 (1)矩阵乘 矩阵乘与数组乘有着较大差别,运算结果也完全不同。矩阵乘的运算符为“*”,运算是按矩阵的乘法规则进行,即参与乘运算的两矩阵的内维必须相同。设A、B为参与乘运算的 =A m×k B k×n。因此,参与运两矩阵,C为A和B的矩阵乘的结果,则它们必须满足关系C m ×n 算的两矩阵的顺序不能任意调换,因为A*B和B*A计算结果很可能是完全不一样的。如:A=[1 1 1;2 2 2;3 3 3]; B=A;

A*B ans= 6 6 6 12 12 12 18 18 18 F=ones(1,3); G=ones(3,1); F*G ans 3 G*F ans= 1 1 1 1 1 1 1 1 1 (2)数组乘 数组乘的运算符为“.*”,运算符中的点号不能遗漏,也不能随意加空格符。参加数组乘运算的两数组的大小必须相等(即同维数组)。数组乘的结果是将两同维数组(矩阵)的对应元素逐一相乘,因此,A.*B和B.*A的计算结果是完全相同的,如: A=[1 1 1 1 1;2 2 2 2 2;3 3 3 3 3]; B=A; A.*B ans= 1 1 1 1 1 4 4 4 4 4 9 9 9 9 9 B.*A ans= 1 1 1 1 1 4 4 4 4 4 9 9 9 9 9 由于矩阵运算和数组运算的差异,能进行数组乘运算的两矩阵,不一定能进行矩阵乘运算。如 A=ones(1,3); B=A; A.*B ans= 1 1 1 A*A ???Error using= =>

数值分析算法在matlab中的实现

数值分析matlab实现高斯消元法: function[RA,RB,n,X]=gaus(A,b) B=[A b];n=length(b);RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0, disp('请注意:因为RA~=RB,所以此方程组无解.') return end if RA==RB if RA==n disp('请注意:因为RA=RB=n,所以此方程组有唯一解.') X=zeros(n,1);C=zeros(1,n+1); for p=1:n-1 for k=p+1:n m=B(k,p)/B(p,p);B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1); end end b=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n); for q=n-1:-1:1 X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q); end else disp('请注意:因为RA=RB0, disp('请注意:因为RA~=RB,所以此方程组无解.') return end if RA==RB if RA==n disp('请注意:因为RA=RB=n,所以此方程组有唯一解.') X=zeros(n,1);C=zeros(1,n+1); for p=1:n-1

数值分析matlab代码

1、%用牛顿法求f(x)=x-sin x 的零点,e=10^(-6) disp('牛顿法'); i=1; n0=180; p0=pi/3; tol=10^(-6); for i=1:n0 p=p0-(p0-sin(p0))/(1-cos(p0)); if abs(p-p0)<=10^(-6) disp('用牛顿法求得方程的根为') disp(p); disp('迭代次数为:') disp(i) break; end p0=p; end if i==n0&&~(abs(p-p0)<=10^(-6)) disp(n0) disp('次牛顿迭代后无法求出方程的解') end 2、disp('Steffensen加速'); p0=pi/3; for i=1:n0 p1=0.5*p0+0.5*cos(p0); p2=0.5*p1+0.5*cos(p1); p=p0-((p1-p0).^2)./(p2-2.*p1+p0); if abs(p-p0)<=10^(-6) disp('用Steffensen加速求得方程的根为') disp(p); disp('迭代次数为:') disp(i) break; end p0=p; end if i==n0&&~(abs(p-p0)<=10^(-6)) disp(n0) disp('次Steffensen加速后无法求出方程的解') end 1、%使用二分法找到方程 600 x^4 -550 x^3 +200 x^2 -20 x -1 =0 在区间[0.1,1]上的根, %误差限为 e=10^-4 disp('二分法')

a=0.2;b=0.26; tol=0.0001; n0=10; fa=600*(a.^4)-550*(a.^3)+200*(a.^2)-20*a-1; for i=1:n0 p=(a+b)/2; fp=600*(p.^4)-550*(p.^3)+200*(p.^2)-20*p-1; if fp==0||(abs((b-a)/2)0 a=p; else b=p; end end if i==n0&&~(fp==0||(abs((b-a)/2)

数值分析幂法与反幂法-matlab程序

数值分析幂法与反幂法 matlab程序 随机产生一对称矩阵,对不同的原点位移和初值(至少取3个)分别使用幂法求计算矩阵的主特征值及主特征向量,用反幂法求计算矩阵的按模最小特征值及特征向量。 要求 1)比较不同的原点位移和初值说明收敛性 2)给出迭代结果,生成DOC文件。 3)程序清单,生成M文件。 解答: >> A=rand(5) %随机产生5*5矩阵求随机矩阵 A = 0.7094 0.1626 0.5853 0.6991 0.1493 0.7547 0.1190 0.2238 0.8909 0.2575 0.2760 0.4984 0.7513 0.9593 0.8407 0.6797 0.9597 0.2551 0.5472 0.2543 0.6551 0.3404 0.5060 0.1386 0.8143 >> B=A+A' %A矩阵和A的转置相加,得到随机对称矩阵B B = 1.4187 0.9173 0.8613 1.3788 0.8044 0.9173 0.2380 0.7222 1.8506 0.5979 0.8613 0.7222 1.5025 1.2144 1.3467 1.3788 1.8506 1.2144 1.0944 0.3929 0.8044 0.5979 1.3467 0.3929 1.6286

B=?? ????? ???? ?? ???6286.13929.03467.15979.08044 .03929.00944 .12144.18506 .13788.13467.12144.15025.17222.08613.05979.08506.17222.02380.09173.08044.03788.18613 .09173 .04187.1 编写幂法、反幂法程序: function [m,u,index,k]=pow(A,u,ep,it_max) % 求矩阵最大特征值的幂法,其中 % A 为矩阵; % ep 为精度要求,缺省为1e-5; % it_max 为最大迭代次数,缺省为100; % m 为绝对值最大的特征值; % u 为对应最大特征值的特征向量; % index ,当index=1时,迭代成功,当index=0时,迭代失败 if nargin<4 it_max=100; end if nargin<3 ep=1e-5; end n=length(A); index=0; k=0; m1=0; m0=0.01; % 修改移位参数,原点移位法加速收敛,为0时,即为幂法 I=eye(n) T=A-m0*I while k<=it_max v=T*u; [vmax,i]=max(abs(v)); m=v(i); u=v/m; if abs(m-m1)

MATLAB与数值分析实验报告一

MATLAB与数值分析实验报告 报告人:秦旸照 学号: 2015020901033 时间: 2016.4.8 电子科技大学电子工程学院

一、实验目的 实验一:MATLAB软件平台与程序设计实验 二、实验原理 1.熟练掌握矩阵的生成、加、减、乘、除、转置、行列式、逆、范数等运算操作。(用.m文件和Matlab函数编写一个对给定矩阵进行运算操作的程序) 2. 熟练掌握算术符号操作和基本运算操作,包括矩阵合并、向量合并、符号转换、展开符号表达式、符号因式分解、符号表达式的化简、代数方程的符号解析解、特征多项式、函数的反函数、函数计算器、微积分、常微分方程的符号解、符号函数的画图等。(用.m文件编写进行符号因式分解和函数求反的程序) 3. 掌握Matlab函数的编写规范。 4.掌握Matlab常用的绘图处理操作,包括:基本平面图、图形注释命令、三维曲线和面的填充、三维等高线等。(用.m文件编写在一个图形窗口上绘制正弦和余弦函数的图形,并给出充分的图形注释) 5. 熟练操作MATLAB软件平台,能利用M文件完成MATLAB的程序设计。 三、实验方案 1. 编程实现以下数列的图像,用户能输入不同的初始值以及系数。并以 x,y为坐标显示图像 x(n+1) = a*x(n)-b*(y(n)-x(n)^2); y(n+1) = b*x(n)+a*(y(n)-x(n)^2) 2. 编程实现奥运5环图,允许用户输入环的直径。 3. 实现对输入任意长度向量元素的冒泡排序的升序排列。 不允许使用sort函数。 四、实验结果 1. 编程实现以下数列的图像,用户能输入不同的初始值以及系数。并以 x,y为坐标显示图像

matlab数值分析例题

1、 在MATLAB 中用Jacobi 迭代法讨论线性方程组, 1231231234748212515 x x x x x x x x x -+=?? -+=-??-++=? (1)给出Jacobi 迭代法的迭代方程,并判定Jacobi 迭代法求解此方程组是否收敛。 (2)若收敛,编程求解该线性方程组。 解(1):A=[4 -1 1;4 -8 1;-2 1 5] %线性方程组系数矩阵 A = 4 -1 1 4 -8 1 -2 1 5 >> D=diag(diag(A)) D = 4 0 0 0 -8 0 0 0 5 >> L=-tril(A,-1) % A 的下三角矩阵 L = 0 0 0 -4 0 0 2 -1 0 >> U=-triu(A,1) % A 的上三角矩阵 U = 0 1 -1 0 0 -1 0 0 0 B=inv(D)*(L+U) % B 为雅可比迭代矩阵 B = 0 0.2500 -0.2500 0.5000 0 0.1250 0.4000 -0.2000 0 >> r=eigs(B,1) %B 的谱半径

r = 0.3347 < 1 Jacobi迭代法收敛。 (2)在matlab上编写程序如下: A=[4 -1 1;4 -8 1;-2 1 5]; >> b=[7 -21 15]'; >> x0=[0 0 0]'; >> [x,k]=jacobi(A,b,x0,1e-7) x = 2.0000 4.0000 3.0000 k = 17 附jacobi迭代法的matlab程序如下: function [x,k]=jacobi(A,b,x0,eps) % 采用Jacobi迭代法求Ax=b的解 % A为系数矩阵 % b为常数向量 % x0为迭代初始向量 % eps为解的精度控制 max1= 300; %默认最多迭代300,超过300次给出警告D=diag(diag(A)); %求A的对角矩阵 L=-tril(A,-1); %求A的下三角阵 U=-triu(A,1); %求A的上三角阵 B=D\(L+U); f=D\b; x=B*x0+f; k=1; %迭代次数 while norm(x-x0)>=eps x0=x; x=B*x0+f; k=k+1; if(k>=max1) disp('迭代超过300次,方程组可能不收敛'); return; end end

数值分析 matlab 实验4

(1) 解题过程如下: (1)MATLAB中创建复化梯形公式和复化辛普森公式的 M 文件:1)复化梯形公式文件: function s=T_fuhua(f,a,b,n) h=(b-a)/n; s=0; for k=1:(n-1) x=a+h*k; s=s+feval(f,x); end s=h*(feval(f,a)+feval(f,b))/2+h*s; 2)复化辛普森公式文件: function s=S_fuhua(f,a,b,n) h=0; h=(b-a)./(2*n); s1=0; https://www.360docs.net/doc/f12321143.html, -5- s2=0; for k=1:n-1 x=a+h*2*k; s1=s1+feval(f,x); end for k=1:n x=a+h*(2*k-1); s2=s2+feval(f,x); end

s=h*(feval(f,a)+feval(f,b)+s1*2+s2*4)/3; 在MATLAB中输入: f=inline('x/(4+x^2)');a=0;b=1; %inline 构造内联函数对象 for n=2:10 s(n-1)=T_fuhua(f,a,b,n);s(n-1)=vpa(s(n-1),10); %调用复化梯形公式,生成任意精度的数值 end exact=int('x/(4+x^2)',0,1);exact=vpa(exact,10) %求出积分的精确值 输出结果:exact = .1115717755 s = Columns 1 through 6 0.1088 0.1104 0.1109 0.1111 0.1113 0.1114 Columns 7 through 9 0.1114 0.1114 0.1115 在MATLAB中输入以下函数用以画出计算误差与 n 之间的曲线: r=abs(exact-s); n=2:10; plot(double(n),double(r(n-1))) 得到结果如图所示: (2)在 MATLAB中输入以下程序代码: f=inline('x/(4+x^2)');a=0;b=1;n=9; %inline 构造内联函数对象 t=T_fuhua(f,a,b,n);t=vpa(t,10) s=S_fuhua(f,a,b,n);s=vpa(s,10)

数值分析实验— MATLAB实现

数值分析实验 ——MATLAB实现 姓名sumnat 学号2013326600000 班级13级应用数学2班 指导老师 2016年1月

一、插值:拉格朗日插值 (1) 1、代码: (1) 2、示例: (1) 二、函数逼近:最佳平方逼近 (2) 1、代码: (2) 2、示例: (2) 三、数值积分:非反常积分的Romberg算法 (3) 1、代码: (3) 2、示例: (4) 四、数值微分:5点法 (5) 1、代码: (5) 2、示例: (6) 五、常微分方程:四阶龙格库塔及Adams加速法 (6) 1、代码:四阶龙格库塔 (6) 2、示例: (7) 3、代码:Adams加速法 (7) 4、示例: (8) 六、方程求根:Aitken 迭代 (8) 1、代码: (8) 2、示例: (9) 七、线性方程组直接法:三角分解 (9) 1、代码: (9) 2、示例: (10) 八、线性方程组迭代法:Jacobi法及G-S法 (11) 1、代码:Jacobi法 (11) 2、示例: (12) 3、代码:G-S法 (12) 4、示例: (12) 九、矩阵的特征值及特征向量:幂法 (13) 1、代码: (13) 2、示例: (13)

一、插值:拉格朗日插值 1、代码: function z=LGIP(x,y)%拉格朗日插值 n=size(x); n=n(2);%计算点的个数 syms a; u=0;%拉格朗日多项式 f=1;%插值基函数 for i=1:n for j=1:n if j==i f=f; else f=f*(a-x(j))/(x(i)-x(j)); end end u=u+y(i)*f;f=1; end z=expand(u);%展开 2、示例: >> x=1:6; y1=x.^5+3*x.^2-6; y2=sin(x)+sqrt(x); >> f1=LGIP(x,y1) f1 = -6+3*a^2+a^5 %可知多项式吻合得很好 >> f2=vpa(LGIP(x,y2),3) f2 = .962e-1*a^4+1.38*a+.300*a^2+.504-.436*a^3-.616e-2*a^5

数理方程基于matlab的数值解法

数理方程数值解法与其在matlab软件上的实现张体强1026222 廖荣发1026226 [摘要] 数学物理方程的数值解在实际生活中越来越使用,首先基于偏微分数值解的思想上,通过matlab软件的功能,研究其数学物理方程的数值解,并通过对精确解和数值解进行对比,追究其数值解的可行性,在此,给出相关例子和程序代码,利于以后的再次研究和直接使用。 [关键字] 偏微分方程数值解matlab 数学物理方程的可视化 一:研究意义 在我们解数学物理方程,理论上求数学物理方程的定解有着多种解法,但是有许多定解问题却不能严格求解,只能用数值方法求出满足实际需要的近似解。而且实际问题往往很复杂,这时即便要解出精确解就很困难,有时甚至不可能,另一方面,在建立数学模型时,我们已作了很多近似,所以求出的精确解也知识推导出的数学问题的精确解,并非真正实际问题的精确解。因此,我们有必要研究近似解法,只要使所求得的近似解与精确解之间的误差在规定的范围内,则仍能满足实际的需要,有限差分法和有限元法是两种最常用的

求解数学物理方程的数值解法,而MATLAB 在这一方面具有超强的数学功能,可以用来求其解。 二:数值解法思想和步骤 2.1:网格剖分 为了用差分方法求解上述问题,将求解区域 {}(,)|01,01x t x t Ω=≤≤≤≤作剖分。将空间区间[0,1]作m 等分,将时 间[0,1]区间作n 等分,并记 1/,1/,,0,,0j k h m n x jh j m t k k n ττ===≤≤=≤≤。分别称h 和τ 为空间和 时间步长。用两簇平行直线,0,,0j k x x j m t t k n =≤≤=≤≤将Ω分割成矩形网格。 2.2:差分格式的建立 0u u t x ??-=??………………………………(1) 设G 是,x t 平面任一有界域,据Green 公式(参考数学物理方程第五章): ( )()G u u dxdt udt udx t x Γ??-=--??? ? 其中G Γ=?。于是可将(1)式写成积分守恒形式: ()0udt udx Γ --=? (2) 我们先从(2)式出发构造熟知的Lax 格式设网格如下图所示

相关文档
最新文档