Matlab软件在数值分析中的应用

《数值分析》课程设计Matlab软件在数值分析中的应用

院（系）名称

2013 年05 月31日

数值分析课程设计评阅书

题目MATLAB软件在数值分析中的应用

学生姓名学号

指导教师评语及成绩

指导教师签名：

年月日答辩评语及成绩

答辩教师签名：

年月日教研室意见

总成绩：

教研室主任签名：

年月日

课程设计任务书

2012—2013学年第二学期

专业班级：学号：姓名：

课程设计名称：数值分析Ⅰ、Ⅱ

设计题目：Matlab软件在数值分析中的应用

完成期限：自2013 年05月21 日至2013年05 月31日共10天

设计依据、要求及主要内容：

一、设计目的

实验方法与理论方法是推进科学技术发展的两大基本方法，由于这些问题不可能用实验手段来实现，因此，学习研究了数值模型.Matlab软件在解决数值模型中起到关键的作用,对于数值分析的研究也十分重要.

二、设计内容

（1）Matlab在数值代数、数值逼近、常微分方程中的应用（2）Matlab在数值分析中的具体应用方法（3）Matlab软件中用不同的方法解决一些数学方面难解的问题．

三、设计要求

1．Matlab软件利用编程解决一些复杂问题，对重点部分加以分析.

2．举一些例子使用所用的方法编写Matlab程序求解.

3．对于在数值分析中的应用应该总结概括完整.

计划答辩时间：2013年06 月 5 日

工作任务与工作量要求：

查阅文献资料不少于3篇，课程设计报告1篇不少于3000字．

指导教师（签字）：教研室主任（签字）：

批准日期：2013 年05 月20 日

Matlab软件在数值分析中的应用

摘要

Matlab已经成为国际上最流行的科学与工程计算的软件工具，是当今最优秀的科学计算软件之一，它集数值计算、图形处理、符号推演、文字处理、动态仿真等功能于一身.Matlab软件在数值分析中的应用也越来越重要,分类讨论了Matlab在数值分析中的应用如：数值代数、数值逼近、常微分方程的计算中的应用,并给出了相应的求解方法.

关键词：数值分析,Matlab,数值代数,数值逼近,微分方程

目录

1 前言 (1)

2 解题思想和方法 (1)

2.1数值方法的特点: (1)

2.2数制与浮点运算 (1)

2.3误差分析 (2)

2.4下面对于范数和极限的讨论 (2)

2.5解线性方程组的直接法 (3)

2.6线性代数方程组的迭代解法 (3)

2.7矩阵特征值问题的求解 (4)

2.8非线性方程求根 (7)

2.9函数插值 (8)

2.11一致逼近与平方逼近 (9)

2.12数值积分 (9)

3 对方程Matlab求解 (10)

3.1幂法的Matlab程序 (10)

3.2 用Newton法 (11)

总结 (13)

参考文献 (13)

1 前言

数值分析能提高自己的思维能力，锻炼自己的理解能力还有学习能力.通过本学期的学习，使我清楚地认识到数值分析对于我们专业的重要性.它也将运用到我们生活中的方方面面，对我们今后的工作生活也有很大的作用.而学习数值分析，我们需要树立算法意识、培养算法设计与分析的能力和应用计算的能力，掌握不同数学问题的求解理论与方法.还有更多的知识需要我们去深入研究，去获取更多的知识.还有就是数值分析在数学建模中也非常重要，他锻炼了一个人的思维能力，在建模中能够较好的去分析问题，更好的去寻找思路.数值分析也离不开Matlab，在学习中我们也常常用Matlab来编程一种计算方式的计算算法.在计算不同数学问题的求解理论方法要树立算法意识和培养算法设计与分析的能力和应用计算机的能力，通过一学期的学习深深感悟到要想把数值分析学习好，需要很深的数学功底和努力学习，在学习中我们要结合实际深刻理解定理、算法，不断有意识和无意识的发现并接受定理，算法蕴含的思想，让它为自我知能的一部分，去引领我们的生活，这样的数学学习才是完美的.

2 解题思想和方法

首先要有可靠的理论分析，以确保算法在理论上的收敛性和数值上的稳定性.

其次要对计算的结果进行误差估计，以确定其是否满足精度，还要考虑算法的运行效率即算法的运算量和存储量.

以及编写Matlab程序实现在计算机上的应用，使其能在一些问题上更快地收敛．2.1数值方法的特点:

1)数值结果要能算得出来，其次结果应有一定的精度

2)误差满足指定的值

3)计算时间应尽可能少

2.2数制与浮点运算

1、数制即在计算机中一个数通常可以用十进制、二进制、八进制等表示.

2、浮点数即用整数部分加小数部分来表示一个实数.

对于二进制的转换举例说明： 10(291)=2(100100011) 2.3误差分析

误差*

x x -的具体数值无法确定，设法给出其绝对值的一个上界

*x x E ε-=≤

绝对误差 *x x ε-=

相对误差 *

x x x

δ-= %x 的相对误差 对于方法的稳定：一种数值方法的传播误差应该可以控制 2.4下面对于范数和极限的讨论

学习它主要目的是利用范数求特征值计算矩阵方程组 1向量范数

11n

j j x x ==∑

{12max ,,,}n x

x x x ∞

=

21/221

()n

j j x x ==∑

2矩阵范数

设.是n R 中的向量范数，对于任何m n R ?，若定义1

max x A Ax ==则

A 是矩阵A 的范数，得到的范数即为向量范数的从属范数.

1范数的从属范数 111

max m

ij j n

i A a <≤==∑

∞范数的从属范数 11

max n

ij i m

j A a ∞<≤==∑

2范数的从属范数 1

(1)

2

n

k n n k =+=

∑ 另：矩阵的F 范数与谱范数均与2范数相容但F 范数不是从属范数，从属范数是所有相容范数中的最小者.

2A ≤ F A

对于任何,0()n n k A R A k ?∈→→∞的充分必要条件是()1A ρ<. 2.5解线性方程组的直接法 1.回迭过程的实现：

1/()/n n nn n

i

i ij j ii j i x b u x b u x u =+=???

=-??

∑ i=n-1,….,1 整个回迭过程乘除法运算量为：

1

(1)

2

n

k n n k =+=

∑. 2.三对角矩阵是一类很重要的特殊矩阵,在数学和物理学中有广泛的应用.三对角矩阵的特征，用待定系数法求解三对角线性方程组的数值解，并与常用的LU 分解法从理论分析和数据实验两方面进行比较,结果表明,两者的时间复杂性前者稍差,而精度两者则相当,最后写出两者的C 程序并运行结果.接下来用一种简单和容易实现的方法求出三对角矩阵的行列式，再利用其逆矩阵可以分解成两个很特殊的矩阵的乘积，给出一种算法实现三对角矩阵的逆的简便计算.

LU 分解：设A 的前n-1个顺序主子阵非奇异，则存在单位下三角矩阵L 及上三角矩阵U 使A=LU 且分解是唯一的. 2.6线性代数方程组的迭代解法

迭代法发的一般形式，Ax b =改写成x Hx g =+，H 为n n R ?矩阵g n R ∈向量而(0)n x R ∈定义向量序列;

(1)()k k x Hx g +=+, {}k x 为迭代序列

解方程组的方法：Jaclbi 迭代法和Gauss-Seidel 迭代法，超松弛迭代法,对于一般的方程组Ax=b,其中n n A R ?∈，detA ≠0,记A=D-L-U.D 为对角部分 L 为严格下三角部分 U 严格上三角矩阵. 设方程组 i x B x f =+

11()j B D L U I D A --=+=- 1J f D b -=

1(1)

()

()11

1 ()i n

k k k i

i ij j ij j j j i ii x b a x a x a -+==+=--∑∑ 此为Jaclbi 迭代法 .

11()()G B D L U I D L A --=-=-- 1()g f D L b -=-

1(1)

(1)()11

1 ()i n

k k k i

i ij j ij j j j i ii x b a x a x a -++==+=--∑∑ i=1,2,…..,n;k=1,2,…. 另外：编程时双保险，避免死循环给出值满足精度，知最大循环次数. 2.7矩阵特征值问题的求解

它包括矩阵特征值和特征向量的计算，而矩阵特征值的计算方法本质上都是迭代法.解决问题的方法有乘幂法和反乘幂法. 乘幂法的基本公式

()(1)

()()()

max()/k k k k k k k

y Az m y z y m -?=?=??=? 带原点位移的反幂法的计算公式

()(1)()

()()

()max()

/k k i k k k k k

A I y z m y z y m λ-?-=?=??=? 在Matlab 中的实验举例说明：

function [k,lambda,vk,wc]=ydwyfmf(A,v0,jlamb,jd,max1) [n,n]=size(A);

A1=A-jlamb*eye(n);jd=jd*0.1;

RAl=det(A1);

if RA1==0

disp('A-aE的n阶行列式n1等于0，所以A-aE不能进行LU分解')

return

end

lambda=0;

if RAl~=0

for p=1:n

h(p)=det(A1(1:p,1:p));

end

hl=h(1:n);

if h(1,i)==0

disp('因为A-aE的各阶主子式等于0，所以A-aE不能进行LU分解')

return

end

end

if h(1,i)~=0

disp('因为A-aE的各阶主子式都不等于0，所以A-aE不能进行LU分解')

k=1;wc=1;state=1;vk=v0;

while((k<=max1)&(state==1))

[L U]=lu(A);

Yk=L/Vk;

Vk=U/Yk;

[m j]=max(abs(Vk));

mk=m;

Vk1=Vk/mk;

Yk1=L/Vk1;

Vk1=U/Yk1;

[m j]=max(abs(Vk1));

mk1=m;

Vk2=(1/mk1)*Vk1;

tzw1=abs((mk-mk1)/mk1);

tzw2=abs(mk1-mk);

Txw1=norm(Vk)-norm(Vk1);

Txw2=(norm(Vk)-norm(Vk1)) rm(Vk1);

Txw=min(Txw1,Txw2);tzw=min(tzw1,tzw2);Vk=Vk2;

mk=mk1;wc=max(Txw,tzw);Vk=Vk2;mk=mk1;state=0;

if(wc>jd)

state=1;

k=k+1;%Vk=Vk2,mk=mk1

end

if(Wc<=jd)

disp('A-aE的秩R（A-aE）和各阶顺序主子式值h1，迭代次数k，按模最小特征值的近似值lambda，特征向量的近似向量Vk，相邻两次迭代的误差Wc如下：') else

disp('A-aE的秩R（A-aE)和各阶顺序主子式值h1，迭代次数k已经达到最大迭代次数max1,按模最小特征值的迭代值lambda，特征向量的迭代向量Vk，相邻两次迭代的误差Wc如下：')

end

h1,RA1

end

end

[V,D]=eig(A,'nobalance');

Vk;k=k-1;Wc;lambdan=jlamb+1/mk1;

2.8非线性方程求根

求()f x =0的解，设*x 为方程的根，即()*0f x =.如果存在正整数使

()*()()m f x x x g x =-，且0< *()g x <∞,则称*x 为他的m 重根.当()f x 在*x 可微时，*

x 为()f x 的重根的充要条件是*x 是()'f x 的根.确定根的存在性定理：若()f x 在[],a b 上连续，且()().0f a f b <，则()f x 在[],a b 上至少有一个根，

steffensen 方法

21

()()k k k k k k k k y f x y x x f x y y +=??

?

=-?+-?

关于它还有一个重要定理，如果()f x 在*x 处二次连续可微，且'*()0f x ≠，则steffensen 方法局部平方收敛，如果'*(1)*()()0m f x f x -=== ，*()0m f x ≠，则s t e f f e n s e n 方法局部线性收敛.

切线法（Newton 法）在Matlab 中编程如下： function [C,D]=newploy(X,Y) n=length(X); D=zeros(n,n); D(:,1)=Y'; for j=2:n

for k=j:n

D(k,j)=(D(k,j-1)-D(k-1,j-1))/(X(k)-X(k-j+1)); end end C=D(n,n); for k=(n-1):-1:1

C=conv(C,poly(X(k)));

m=length(C); C(m)=C(m)+D(k,k); end 2.9函数插值

插值问题的提法是：假定区间[a ，b]上的实值函数f （x ）在该区间上 n+1个互不相同点x 0，x 1……x n 处的值是f (x 0)，……f（x n ），要求估算f （x ）在[a ，b]中某点x*的值.基本思路是，找到一个函数P(x)，在x 0，x 1……x n 的节点上与f(x)函数值相同(有时，甚至一阶导数值也相同)，用P(x*)的值作为函数f(x*)的近似.拉格朗日插值公式对给定的n+1个节点x 0,x 1,x 2,…,x n 及对应的函数值y 0,y 1,y 2,…,y n ， 构造一个n 次插值多项式：

)(0x l y y k n

k k ∑==

简单程序示例:

x=[0.0 0.1 0.195 0.3 0.401 0.5];

y=[0.39849 0.39695 0.39142 0.38138 0.368120.35206 ]; T=interp1(x,y,0..25,'linear') %线性插值 (返回结果T=0.3862)

T=interp1(x,y,0.25,'nearest') % 两点插值 (返回结果T=0.3814)

T=interp1(x,y,0.25,'spline') % 三次样条插值 (返回结果T =0.3867)

T=interp1(x,y,0.25,'cubic') %三次插值 (返回结果T =0.3867) 2.10样条函数

在插值问题中，样条插值通常比多项式插值好用.用低阶的样条插值能产生和高阶的多项式插值类似的效果，并且可以避免被称为龙格现象的数值不稳定的出现.并且低阶的样条插值还具有“保凸”的重要性质.

2.11一致逼近与平方逼近

近似代替又称为逼近，函数 f (x ) 称为被逼近函数；P (x ) 称为逼近函数，两者之差称为逼近的误差.

1、一致逼近（均匀逼近）

以)()(max x p x f b

x a -≤≤ 作为度量误差f (x )- P (x ) 的“大小” 标准.

2、 平方逼近（均方逼近）

以[]12

2

(()())b

a f x p x dx -?作为度量误差f (x )- P (x )的“大小” 标准.

function c=pf1(f,n,a,b) %c 最佳平方逼近多项式的系数

syms x; for i=1:n+1 for j=1:n+1

A(i,j)=int(x^(i+j-2),x,a,b); end

B(i)=int(f*x^(i-1),x,a,b); end c=inv(A)*B';

2.12数值积分

1.插值型的求积公式:构造出的求积0

()n

n k k k f x =I =A ∑,()b

k k a

l x dx A =? 公式积称为是插值

型求积公式

2.牛顿——柯特斯公式：设将积分区间划分为n 等分，步长 b a

h n

-=

，选取等距节点 k x a kh =+ ,构造出的插值型求积公式:()0

()()n

n n k k k b a C f x =I =-∑

称为牛顿—柯特斯公式，式中

()001n

n n k

j j k

t j C

dt n k j

=≠-=-∏?

例

dx x x

?+1

024

x=0:1/8:1; y=x./(4+x.^2); trapz(x,y) 运行结果：

ans= 0.11140235452955 3 利用Matlab 求解

3.1幂法的Matlab 程序

用幂法计算矩阵的主特征值和对应的特征向量的Matlab 主程序 function [k,lambda,Vk,Wc]=mifa(A,V0,jd,max1) lambda=0;k=1;Wc =1; ,jd=jd*0.1;state=1; V=V0; while((k<=max1)&(state==1)) Vk=A*V; [m j]=max(abs(Vk)); mk=m; tzw=abs(lambda-mk); Vk=(1/mk)*Vk; Txw=norm(V-Vk); Wc=max(Txw,tzw); V=Vk;lambda=mk;state=0; if(Wc>jd) state=1; end

k=k+1;Wc=Wc; end if(Wc<=jd)

disp('请注意：迭代次数k,主特征值的近似值lambda,主特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下：')

else

disp('请注意：迭代次数k 已经达到最大迭代次数max1,主特征值的迭代值

lambda,主特征向量的迭代向量Vk,相邻两次迭代的误差Wc 如下：')

end

Vk=V;k=k-1;Wc;

用幂法计算下列矩阵的主特征值和对应的特征向量的近似向量，精度. 解 （1）输入MATLAB 程序 A=[1 -1;2 4];

V0=[1,1]'；

[k,lambda,Vk,Wc]=mifa(A,V0,0.00001,100),

[V ,D] = eig (A), Dzd=max(diag(D)), wuD= abs(Dzd- lambda), wuV=V(:,2)./Vk, 运行后屏幕显示结果

请注意：迭代次数k,主特征值的近似值lambda,主特征向量的近似向量Vk,相邻两次迭代的误差Wc 如下： k = lambda = Wc =

33 3.00000173836804 8.691862856124999e-007

Vk = V = wuV =

-0.49999942054432

-0.70710678118655

0.44721359549996

-0.89442822756294

1.00000000000000

0.70710678118655

-0.89442719099992

-0.89442719099992

Dzd = wuD =

3 1.738368038406435e-006

由输出结果可看出，迭代33次，相邻两次迭代的误差Wc 8.69 19e-007，矩阵的主特征值的近似值lambda3.000 00和对应的特征向量的近似向量Vk （-0.500 00，1.000 00, lambda 与例5.1.1中的最大特征值近似相等，绝对误差约为1.738 37e-006，Vk 与特征向量 的第1个分量的绝对误差约等于0，第2个分量的绝对值相同.由wuV 可以看出，的特征向量V(:,2) 与Vk 的对应分量的比值近似相等.因此，用程序mifa.m 计算的结果达到预先给定的精度. 3.2 用Newton 法

例.使用Newton 法计算方程0107222=--x e x 在区间[]2,2-上的根，其初始值

x．

5.1

解：编写主程序为：

[x,time]=Newton('2*exp(2*x)-7*x^2-10',1.5,1e-5)

x=0:0.01:2;

y=2*exp(2*x)-7*x.^2-10;

plot(x,y)

grid

运行并调用Newton的主程序结果为：

wucha=0.2430

wucha=0.0615

wucha=0.0311

……………

wucha=2.3876e-005

wucha=1.5073e-005

wucha=9.5154e-006

x=1.1197

time=20

总结

通过一周的时间，数值分析课程设计终于完成了．在这期间，我收获了很多，也学习很多．以前仅仅只学习书本上的理论知识，从没有自己的动手实践.通过做这次的课程设计，使我把理论的知识与实践结合到了一起，达到了理论与实践的结合．此外，也提高了自己的动手能力．而最大难题就是程序，由于不擅长编程序，在编写程序是总是出现问题，于是，自己在网上搜了了一些资料，并通过查阅此资料．并且通过自己的修改终于能运行出来了时候，心里有一点小小的自豪感．

一个学期的数值分析学习，使我有很多感悟.其中我感受最深的是数值分析是一门重视算法和原理的学科，它的内容更接近于实际，像求解线性方程组、矩阵的特征值等，是数学理论更加有实际意义，随着学习的深入，我发现要想学会数值分析必须把那些理论和定理搞得非常透彻，虽然我现在还无法真正理解那些思想，但我会不断的有意识的无意识的去熟悉那些定理和蕴含的思想.进而真正让数学帮助我们全方位的成长.

通过这次课程设计，我从中学会了很多，也发现自己真的还有很多不足以及很多东西需要去学习．所以在以后的生活学习中要不断的扩大自己的视野，多学习一些与专业有关的知识，不能只满足于课本上的知识．所以在完成本专业的基础上，要不断涉猎，完善自我，希望自己在以后的课程中会得到更好的锻练．总的来说这次课程设计还是有很多的收获的，并且特别感谢我们组的成员在做课程设计的过程中对我的帮助．

参考文献

[1] 冯国忱黄明游．数值分析(上册)[M]．北京：高等教育出版社，2007．

[2] 宋叶志贾东志．MATLAB数值分析与应用[M]．北京：机械工业出版社，2009．

[3] 张德丰．MATLAB数值分析与应用[M]．北京：国防工业出版社，2007．

[4] 陈东彦李冬梅，王树忠.数学建模[M].北京：科学出版社，2007

[5]李庆扬王能超易大义.数值分析[M]．4版.北京清华大学出版社,2005

《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实现

第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与数值分析作业

数值分析作业(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．求数值积分： fx=@(x)exp(1./x); >> quadl(fx,1,5) 2．获取x=xlsread('oillack.xls','sheet1','a1:a73') excel文件名是oillack.xls，sheet1是表名，a1:a73'是a列的1到73行 long x=xlsread('F:\学习\大三\大三下\巷道力学模型\新建文件夹(2)\1.xlsx','sheet1','a2:a') 3. 在matlab的图中插入文本框后将文本框旋转的方法： text(0.5,0.6,'渗透率/mD','Rotation',90) 4. matlab中插入一条直线的方法： line([0.01 0.01],[0 1.75]) 5．Matlab 中画三维图 x=-7.5:0.5:7.5; y=x; % 先产生x及y二个阵列 >> [x,y]=meshgrid(x,y); % 再以meshgrid形成二维的网格数据 >> z=x.^2+y.^2; % 产生z轴的数据 >> mesh(x,y,z) % 将z轴的变化值以网格方式画出 >> surf(x,y,z) % 将z轴的变化值以曲面方式画出 Matlab指数拟合方法 x=[1982 1992 2002]； y=[103.5 34.5 23.3]； cftool（x，y） 在弹出的对话框选择fitting，弹出新的对话框选择new fit，然后在第三个下拉菜单（Type of fit）中选择Exponential，然后点击Apply，即可；最后结果 General model Exp1: f(x) = a*exp(b*x) Coefficients (with 95% confidence bounds): a = 1.453e+082 (-7.288e+084, 7.317e+084) b = -0.09312 (-0.3464, 0.1602)

数值分析 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.wendangku.net/doc/de465684.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 格式设网格如下图所示

数值分析五个题目的C语言及Matlab程序

本文档包含上一个文档中的五个数值分析实验题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};

数值分析在生活中的应用举例及Matlab实现

Matlab 实验报告 学院：数学与信息科学学院班级：信息班 学号：20135034027 姓名：马永杉

最小二乘法，用MATLAB实现 1.数值实例 下面给定的是郑州最近1个月早晨7：00左右的天气预报所得到的温度，按照数据找出任意次曲线拟合方程和它的图像。下面用MATLAB编程对上述数据进行最小二乘拟合。 2、程序代码 x=[1:1:30]; y=[9,10,11,12,13,14,13,12,11,9,10,11,12,13,14,12,11,10,9,8,7,8,9,11,9 ,7,6,5,3,1]; a1=polyfit(x,y,3) %三次多项式拟合% a2= polyfit(x,y,9) %九次多项式拟合% a3= polyfit(x,y,15) %十五次多项式拟合% b1=polyval(a1,x) b2=polyval(a2,x) b3=polyval(a3,x) r1= sum((y-b1).^2) %三次多项式误差平方和% r2= sum((y-b2).^2) %九次次多项式误差平方和% r3= sum((y-b3).^2) %十五次多项式误差平方和% plot(x,y,'*') %用*画出x,y图像% hold on plot(x,b1, 'r') %用红色线画出x,b1图像% hold on plot(x,b2, 'g') %用绿色线画出x,b2图像% hold on plot(x,b3, 'b:o') %用蓝色o线画出x,b3图像% 2.流程图

4．数值结果分析 不同次数多项式拟合误差平方和为： r1=67.6659 r2=20.1060 r3=3.7952 r1、r2、r3分别表示三次、九次、十五次多项式误差平方和。 5、拟合曲线如下图