MATLAB应用实例分析例分析

MATLAB应用实例分析例分析

Matlab应用例题选讲

仅举一些运用MATLAB的例子，这些问题在数学建模中时常遇到，希望能帮助同学们在短时间内方

便、快捷的使用MATLAB 解决数学建模中的问题，并善用这一工具。常用控制命令:

clc:%清屏; clear:%清变量; save:%保存变量; load:%导入变量

一、利用公式直接进行赋值计算

本金P以每年n次，每次i%的增值率(n与i的乘积为每年增值额的百分比)增加，当增加到r×P 时所花费的时间T为:(利用复利计息公式可得到下式) lnrnT() r，P,P(1，0.01i),T,r,2,i,0.5,n,12nln(1，0.01i)

MATLAB 的表达形式及结果如下:

>> r=2;i=0.5;n=12; %变量赋值

>> T=log(r)/(n*log(1+0.01*i)) 计算结果显示为:

T = 11.5813

即所花费的时间为T=11.5813 年。

分析:上面的问题是一个利用公式直接进行赋值计算问题，实际中若变量在某个范围变化取很多值时，

使用MATLAB，将倍感方便，轻松得到结果，其绘图功能还能将结果轻松的显示出来，变量之间的变化

规律将一目了然。

若r在[1,9]变化，i在[0.5,3.5]变化;我们将MATLAB的表达式作如下改动，结果如图1。 r=1:0.5:9;

i=0.5:0.5:3.5;

n=12;

p=1./(n*log(1+0.01*i));

T=log(r')*p;

plot(r,T)

xlabel('r') %给x轴加标题

ylabel('T') %给y轴加标题

q=ones(1,length(i));

text(7*q-0.2,[T(14,1:5)+0.5,T(14,6)-0.1,T(14,7)-0.9],num2str(i'))

40

350.530

25

20T 115

1.510 2

2.55 3

3.50123456789r

图1

1

从图1中既可以看到T随r的变化规律，而且还能看到i的不同取值对T—r 曲线的影响(图中的六条曲线分别代表i的不同取值)。

二、已知多项式求根

65432已知多项式为，求其根。 h,x,10x，31x,10x,116x，200x,96

分析:对多项式求根问题，我们常用roots()函数。MATLAB 的表达形式及结果如下: >> h=roots([1 -10 31 -10 -116 200 -96]) %中括号内为多项式系数由高阶到常数。计算结果显示为(其中i为虚数单位):

h =

-2.0000

4.0000

3.0000

2.0000 + 0.0000i

2.0000 - 0.0000i

1.0000

如果已知多项式的根，求多项式，用poly()函数。对上面得到的h的值求多项式，其MATLAB的表达形式及结果如下:

>>h=[-2.0000 4.0000 3.0000 2.0000+0.0000i 2.0000-0.0000i 1.0000];

>>c=poly(h)

计算结果显示为:

c =

1 -10 31 -10 -116 200 -96

三、方程组的求解

8x，x，6x,7.5,123,求解下面的方程组: 3x，5x，7x,4,123,4x，9x，

2x,12123,

分析:对于线性方程组求解，常用线性代数的方法，把方程组转化为矩阵进行计算。

,1 ,x,abax,b,x,a\b

MATLAB 的表达形式及结果如下:

>> a=[8 1 6;3 5 7;4 9 2]; %建立系数矩阵

>> b=[7.5;4;12]; %建立常数项矩阵

>> x=a\b %求方程组的解

计算结果显示为:

x =

1.2931

0.8972

-0.6236

四、数据拟合与二维绘图

在数学建模竞赛中，我们常会遇到这种数据表格问题，如果我们仅凭眼睛观察，很难看到其中的规律，也就更难写出有效的数学表达式从而建立数学模型。因此可以利用MATLAB的拟合函数，即polyfit()

函数，并结合MATLAB的绘图功能(利用plot()函数)，得到直观的表示。

2

例:在化学反应中，为研究某化合物的浓度随时间的变化规律，测得一组数据如下表:

1 2 3 4 5 6 7 8 T(分)

y 4 6.4 8.0 8.4 9.28 9.5 9.7 9.86

9 10 11 12 13 14 15 16 T(分)

y 10 10.2 10.32 10.42 10.5 10.55 10.58 10.6

分析:

MATLAB 的表达形式如下:

t=[1:16]; %数据输入

y=[4 6.4 8 8.4 9.28 9.5 9.7 9.86 10 10.2 10.32 10.42 10.5 10.55

10.58 10.6];

plot(t,y,'o') %画散点图

p=polyfit(t,y,2) %二次多项式拟合

hold on

xi=linspace(0,16,160); %在[0,16]等间距取160 个点

yi=polyval(p,xi); %由拟合得到的多项式及xi，确定yi

plot(xi,yi) %画拟合曲线图

执行程序得到图2;

11

10

9

8

7

6

5

40246810121416

图2

显示的结果为

p=

0.0445 1.0711 4.3252 -

2p的值表示二阶拟合得到的多项式为:y= -0.0445t+1.0711t+ 4.3252 下面是用lsqcurvefit()函数，即最小二乘拟合方法的Matlab表达: t=[1:16];

y=[4 6.4 8 8.4 9.28 9.5 9.7 9.86 10 10.2 10.32 10.42 10.5 10.55 10.58 10.6];

x0=[0.1,0.1,0.1];

zuixiao=inline('x(1)*t.^2+x(2)*t+x(3)','x','t');

x=lsqcurvefit(zuixiao,x0,t,y) %利用最小二乘拟合

其显示的结果为:

x =

-0.0445 1.0711 4.3252

可以看出其得到的结果与polyfit函数的结果相同。这说明在多项式拟合问题

上这两个函数的效果是相同的。

3

下面的一个例子将体现lsqcurvefit()函数的优势。

例2: 在物理学中，为研究某种材料应力与应变的关系，测得一组数据如下表: 925 1125 1625 2125 2625 3125 3625 应力σ

0.11 0.16 0.35 0.48 0.61 0.71 0.85 应变ε

如果假定应力与应变有如下关系(σ为应力值，ε为应变值):ε=a+blnσ试

计算a 、b 的值。

MATLAB 的表达形式如下:

x=[925,1125,1625,2125,2625,3125,3625];

y=[0.11,0.16,0.35,0.48,0.61,0.71,0.85]; plot(x,y,'o')

[p,resid1]=polyfit(x,y,2)

hold on

xi=linspace(700,3700,3000); yi=polyval(p,xi);

plot(xi,yi)

x0=[0.1,0.1];

fff=inline('a(1)+a(2)*log(x)','a','x');

[a,resid2]=lsqcurvefit(fff,x0,x,y) plot(xi,fff(a,xi),'r') 执行程序得到图3，图中蓝色曲线为利用polyfit()函数得到的曲线，红色曲

线为利用lsqcurvefit()函

数得到的曲线;

0.90.80.70.60.50.40.30.20.10-0.15001000150020002500300035004000

其显示的结果为:

p =

-0.0000 0.0004 -0.2266 resid1 =

R: [3x3 double]

df: 4

normr: 0.0331

a =

-3.5810 0.5344

resid2 =

0.0064

其中a的值代表利用lsqcurvefit()函数得到的关系为:ε=-3.5810+0.5344σ

4

resid1、resid2 分别代表运用polyfit()函数、lsqcurvefit()函数得到的残差。可以看出利用lsqcurvefit()函数残差更小，即得到了更好的拟合效果。

在数学建模的实际问题中，如果问题的机理不明，我们只能采用polyfit()函数，即多项式拟合的

方法，以获得近似的数据描述函数;但如果通过分析，可以得到一些机理，那么采用最小二乘的方法将

得到更好的效果，而且得到的拟合函数也更有意义。

五、隐函数的图形绘制

1plot()只能绘制显函数图形，对于形如的复杂隐函数，很难转化为,lny，

ln(,1，y)，x,sin(x),0y

显函数并利用plot()函数绘制图形，这时就可以用ezplot()函数直接绘制其

曲线。 MATLAB的表达形式如下:

>> ezplot('1/y-log(y)+log(-1+y)+x-sin(x)')

执行程序得到图5

x = sin(3 t) cos(t), y = sin(3 t) sin(t)1/y-log(y)+log(-1+y)+x-sin(x) = 060.4

40.2

20

-0.2yy0

-0.4-2-0.6-4-0.8

-6-0.8-0.6-0.4-0.200.20.40.60.8-6-4-20246xx

图5 图6 如果是形如下面的参数方程,同样可以利用ezplot()函数绘制其曲

x,sin3tcost,y,sin3tsint,t,(0,,)线。MATLAB的表达形式如下:

>> ezplot('sin(3*t)*cos(t)','sin(3*t)*sin(t)',[0,pi])

执行程序得到图6。

六、三维图形绘制

假设有一个时间向量t，对该向量进行下列运算则可以构成三个坐标值向量

x,sint,y,cost,z,t

对于上面的方程可以利用ezplot3()函数或plot3()函数绘制三维曲线。这里

仅列举ezplot3()函数的使用。

MATLAB 的表达形式如下:

>> ezplot3('sin(t)','cos(t)','t',[0,6*pi])

执行程序得到图7:

3绘制下述曲面: z(r,,),rcos(3,),其中0,r,1,0,,,2,

MATLAB 的表达形式如下:

nr=12;nth=50;

r=linspace(0,1,nr);

theta=linspace(0,2*pi,nth);

[R,T]=meshgrid(r,theta)

x=cos(theta')*r;

5

y=sin(theta')*r;

surf(x,y,R.^3.*cos(3*T))

执行程序得到图8。

x = sin(t), y = cos(t), z = t

20

15

10z

5

010.510.500-0.5-0.5-1-1yx

图7 图8

除了surf()函数还有surfc()、surfl()、mesh()、waterfall()函数也用于曲面的绘制，具体效果如图9所示，可

以针对自己的需要选取适合的曲面绘制函数。

图9

MATLAB作图之三维绘图示例山体绘制

%三维绘图示例山体绘制

%mesh函数演示

x=1.0:0.1:2.0;

y=2.0:0.1:3.0;

[X,Y]=meshgrid(x,y);

z=[5.11 5.13 5.14 5.13 5.09 5.04 4.98 4.93 4.89 4.85 4.85

5.39 5.49 5.51 5.46 5.32 5.14 4.94 4.74 4.59 4.49 4.48

5.61 5.77 5.81 5.71 5.51 5.23 4.90 4.59 4.36 4.21 4.19

5.73 5.92 5.97 5.86 5.62 5.27 4.88 4.51 4.23 4.05 4.03

5.74 5.92 5.97 5.86 7.62 5.27 4.88 4.51 4.21 4.04 4.02

5.63 5.79 5.84

6.74 10.53 9.23 8.91 4.59 4.33 4.18 4.16

5.42 5.53 5.56 5.49 7.35 5.16 4.93 4.73 4.55 4.45 4.44

5.14 5.18 5.19 5.17 11.12 5.05 4.97 4.90 4.84 4.81 4.80

4.48 4.80 4.79 4.82 4.87 4.94

5.02 5.10 5.16 5.19 5.20

4.56 4.45 4.43 4.49 4.64 4.84

5.06 5.28 5.45 5.55 5.56

4.36 4.19 4.16 4.25 4.47 4.76

5.09 5.41 5.66 5.81 5.83]; figure(1)

%画网格图

mesh(X,Y,z);

6

colormap([0 1 0]);

xlabel('x轴');

ylabel('y轴');

zlabel('z轴');

%画表面图

figure(2)

surf(X,Y,z);

colormap([1 0 0]);

xlabel('x轴');

ylabel('y轴');

zlabel('z轴');

1212

1010

88z轴z轴66

4433221.81.82.52.51.61.61.41.41.21.22211y轴y轴x轴x轴

七、二项分布的使用

飞机成功起飞的概率问题:

由16 架飞机组成的空军飞行中队要求做好立即起飞的准备，其中一架飞机不能立即起飞的概率为

20%，重新起飞需几分钟的时间，因此一架飞机立刻起飞的概率为0.80。12 架飞机能够成功起飞的概率

为多少,

分析:这是一个概率中的二项分布问题，常用binopdf()函数。

h=binopdf(12,16,0.80) %二项分布函数的概率值

计算结果为

h=0.2001

另一方面，至少有14架飞机立刻成功起飞的概率为: h=1-

binocdf(13,16,0.80) % 或 h=sum(binopdf(14:16,16,0.80))，其中binocdf()为二项分布的累积概率值，计算结果为h=0.3518。

在实际的数学建模竞赛中，仅罗列一个一个的数据是枯燥而又不直观的，很难吸引人们的注意，也不容

易打动评委们的心;因此，结合数值计算结果，并合理的利用MATLAB的绘图功能会起到事半功倍的效

果。下面的程序为运行结果的绘图(图10)表示:

n=1:16;

h=binopdf(n,16,0.80);

plot([n;n],[zeros(1,16);h],'k') %二维绘图函数

text(8-.7:16-.7,h(8:16)+.005,num2str(h(8:16)',3)) %在图中进行注释函数 axis([0 17 0 0.27]) %坐标轴取值范围函数

xlabel('Number of aircraft launched on time') %给x 轴加标题

ylabel('probability') %给y 轴加标题

set(gca,'XTick',0:2:16)

7

set(gca,'XTickLabel',{'0 架','2 架','4 架','6 架','8 架','10 架','12 架','14 架','16 架'})

13

0.2460.2512

11 0.211 0.20.210

90.158 0.12 0.113probability70.16

0.05550.05 0.02814 0.01970.005533002468101214160 架2 架4 架6 架8 架10 架12 架14 架16 架Number of aircraft launched on time 图10 图11 根据例5中的数据，用线性回归的方法画出95%的置信区间，见图11，程序如下: x=[1:16]; %数据输入

y=[4 6.4 8 8.4 9.28 9.5 9.7 9.86 10 10.2 10.32 10.42 10.5 10.55 10.58 10.6];

[p,resid1]=polyfit(x,y,2) ;

[yhat,w]=polyconf(p,x,resid1,0.05); plot(x,yhat,'k-',x,yhat-w,'k--',x,yhat+w,'k--',x,y,'ks',[x;x],[yhat;y],'k-')

残差的进一步分析:

首先计算残差值，然后用函数normplot 将其画在图12中，以观察残差是否服从正态分布。程序为: x=[1:16]; %数据输入

y=[4 6.4 8 8.4 9.28 9.5 9.7 9.86 10 10.2 10.32 10.42 10.5 10.55 10.58 10.6];

[p,resid1]=polyfit(x,y,2)

normplot(y-polyval(p,x))

whitebg('white') %设置背景为白色

90

Normal Probability Plot800.98

700.95

0.90 60

0.75 50

400.50 Probability300.25

200.10 100.05

0.02 0155160165170175180185190195-1-0.500.5Data

图12 图13 由图12看出，这些残差点与直线非常接近，由此得出结论，这些残差值非常接近于正态分布，选择的模型是合适的。

8

第06章_MATLAB数值计算_例题源程序汇总

第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;

同济大学数值分析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

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

第3章 MATLAB数值计算-习题 答案

roots([1 -1 -1]) x=linspace(0,2*pi,10); y=sin(x); xi=linspace(0,2*pi,100); y1=interp1(x,y,xi); y2=interp1(x,y,xi,'spline'); y3=interp1(x,y,xi,'cublic'); plot(x,y,'o',xi,y1,xi,y2,xi,y3) x=[0 300 600 1000 1500 2000]; y=[0.9689 0.9322 0.8969 0.8519 0.7989 0.7491]; xi=linspace(0,2000,20); yi=1.0332*exp(-(xi+500)/7756); y1=interp1(x,y,xi,'spline'); subplot(2,1,1);plot(x,y,'o',xi,yi,xi,y1,'*') p=polyfit(x,y,2); y2=polyval(p,xi); subplot(2,1,2);plot(x,y,'o',xi,yi,xi,y2,'*') x=[0 300 600 1000 1500 2000]; y=[0.9689 0.9322 0.8969 0.8519 0.7989 0.7491]; xi=linspace(0,2000,20); y1=interp1(x,y,xi,'spline'); subplot(2,1,1);plot(x,y,'-o', xi,y1,'-*') p=polyfit(x,y,2); y2=polyval(p,xi); subplot(2,1,2);plot(x,y,'-o',xi,y2,'-*')

东南大学-数值分析上机题作业-MATLAB版

2015.1.9 上机作业题报告 JONMMX 2000

1.Chapter 1 1.1题目 设S N =∑1j 2?1 N j=2 ，其精确值为 )1 1 123(21+--N N 。 （1）编制按从大到小的顺序1 1 131121222-+ ??+-+-=N S N ，计算S N 的通用程序。 （2）编制按从小到大的顺序1 21 1)1(111222-+ ??+--+-= N N S N ，计算S N 的通用程序。 （3）按两种顺序分别计算64210,10,10S S S ,并指出有效位数。（编制程序时用单精度） （4）通过本次上机题，你明白了什么？ 1.2程序 1.3运行结果

1.4结果分析 按从大到小的顺序，有效位数分别为：6，4，3。 按从小到大的顺序，有效位数分别为：5，6，6。 可以看出，不同的算法造成的误差限是不同的，好的算法可以让结果更加精确。当采用从大到小的顺序累加的算法时，误差限随着N 的增大而增大，可见在累加的过程中，误差在放大，造成结果的误差较大。因此，采取从小到大的顺序累加得到的结果更加精确。 2.Chapter 2 2.1题目 （1）给定初值0x 及容许误差ε，编制牛顿法解方程f(x)=0的通用程序。 （2）给定方程03 )(3 =-=x x x f ,易知其有三个根3,0,3321= *=*-=*x x x ○1由牛顿方法的局部收敛性可知存在,0>δ当),(0δδ+-∈x 时，Newton 迭代序列收敛于根x2*。试确定尽可能大的δ。 ○2试取若干初始值，观察当),1(),1,(),,(),,1(),1,(0+∞+-----∞∈δδδδx 时Newton 序列的收敛性以及收敛于哪一个根。 （3）通过本上机题，你明白了什么？ 2.2程序

Matlab大数值计算题目

Matlab大数值计算题目 1、统计附件1中的数据，对其中的数据划分区间，从0到50，每 10个单位一个区间，分为5个区间，统计每个区间的数量，画出柱状图。 Matlab程序： clear;clc;close all Data=xlsread('数据.xls'); Q=0:10:50; n=length(Data); m=length(Q); T=zeros(size(Q)); for s=1:n for t=1:m-1 if Data(s)>Q(t)&Data(s)

2、统计附件2中第二列数据中1至100每个数字出现的总次数， 附件2中第三列为每出现第二列数字所对应的次数，最后画出柱状图。 Matlab程序： clear;clc;close all Data=load('WEIBOIDWITHCOMMENTS.txt'); DATA=Data(:,2); t=Data(:,3); % m=max(DATA); m=100; T=zeros(m,1); for i=1:m data=DATA; data(data~=ones(size(data))*i)=0; data(data~=0)=1; n=data.*t; N=sum(n); T(i)=N; end bar(T)

3、 找到矩阵迷宫的通路,矩阵第1行第1列为迷宫的入口，第8行 第8列为迷宫的出口。（0表示路，1表示墙） 000 00000011 11010000 01010010 11010010 11010010 00011010 01000011 11110?????????????????????????? Matlab 程序： 主程序： clear all clc maze=[0,0,0,0,0,0,0,0; 0,1,1,1,1,0,1,0; 0,0,0,0,1,0,1,0; 0,1,0,0,0,0,1,0; 0,1,0,1,1,0,1,0; 0,1,0,0,0,0,1,1; 0,1,0,0,1,0,0,0; 0,1,1,1,1,1,1,0]; total=0; maze(1,1)=3;

数值分析上机题(matlab版)(东南大学)

数值分析上机题(matlab版)(东南大学)

数值分析上机报告

第一章 一、题目 精确值为)1 1 123(21+--N N 。 1) 编制按从大到小的顺序 1 1 131121222-+??+-+-= N S N ，计算S N 的通用程序。 2) 编制按从小到大的顺序 1 21 1)1(111222-+??+--+-= N N S N ，计算S N 的通用程序。 3) 按两种顺序分别计算6 42 10,10, 10S S S ,并指出有效位 数。（编制程序时用单精度） 4) 通过本次上机题，你明白了什么？ 二、通用程序 clear N=input('Please Input an N (N>1):'); AccurateValue=single((0-1/(N+1)-1/N+3/2)/2); Sn1=single(0); for a=2:N; Sn1=Sn1+1/(a^2-1); end Sn2=single(0); for a=2:N; Sn2=Sn2+1/((N-a+2)^2-1); end fprintf('The value of Sn using different algorithms (N=%d)\n',N); disp('____________________________________________________') fprintf('Accurate Calculation %f\n',AccurateValue); fprintf('Caculate from large to small %f\n',Sn1); fprintf('Caculate from small to large %f\n',Sn2);

MATLAB解决数值分析问题

1. 使用444（x）对数据进行插值，并写出误差分析理论。 建立脚本 x1=[0.2 0.4 0.6 0.8 1.0]; y1=[0.98 0.92 0.81 0.64 0.38]; n=length(y1); c=y1(:); for j=2:n %求差商 for i=n:-1:j c(i)=(c(i)-c(i-1))/(x1(i)-x1(i-j+1)); end end syms x df d; df(1)=1;d(1)=y1(1); for i=2:n %求牛顿差值多项式 df(i)=df(i-1)*(x-x1(i-1)); d(i)=c(i)*df(i); end disp('4次牛顿插值多项式'); P4=vpa(collect((sum(d))),5) %P4即为4次牛顿插值多项式,保留小数点后5位数figure ezplot(P4,[0.2,1.08]); 输出结果为 P4 =- 0.52083*x^4 + 0.83333*x^3 - 1.1042*x^2 + 0.19167*x + 0.98 插值余项：R4(x)=f(5)( ξ)/ (5!)* (x - 0.6)*(x - 0.4)*(x - 0.8)*(x - 1)*(x-0.2) 新建一个M-file

syms x l; x1=[0.2 0.4 0.6 0.8 1.0]; y1=[0.98 0.92 0.81 0.64 0.38]; n=length(x1); Ls=sym(0); for i=1:n l=sym(y1(i)); for k=1:i-1 l=l*(x-x1(k))/(x1(i)-x1(k)); end for k=i+1:n l=l*(x-x1(k))/(x1(i)-x1(k)); end Ls=Ls+l; end Ls=simplify(Ls) %为所求插值多项式Ls(x). figure ezplot(Ls,[0.2,1.08]); 输出结果为 Ls =- (25*x^4)/48 + (5*x^3)/6 - (53*x^2)/48 + (23*x)/120 + 49/50 插值余项：R4(x)=f(5)( ξ)/ (5!)* (x - 0.6)*(x - 0.4)*(x - 0.8)*(x - 1)*(x-0.2) 2. 试求3次、 合曲线，用图示数据曲线及相应的三种拟合曲线。 建立脚本 X=[0.0 0.1 0.2 0.3 0.5 0.8 1.0]; Y=[1.0 0.41 0.50 0.61 0.91 2.02 2.46]; p1=polyfit(X,Y,3) p2=polyfit(X,Y,4)

matlab程序题

1.采用数值计算方法，画出dt t t x y x ?=0sin )(在]10 ,0[区间曲线，并计算)5.4(y 。 （提示：cumtrapz 快捷，在精度要求不高处可用；quad 也可试。巧用find 。）第一种答： 程序代码 clear; t=0:0.1:10; y=sin(t)./t; s =cumtrapz(t,y); plot(t,y) > y='sin(t)./t'; quad(y,0,4.5,0.01) 运行结果： ans = 1.6541 第二种答： 程序代码 d=0.5; %分区间的区间间隔 tt=0:d:10; t=tt+(tt==0)*eps; %使0被一个极小数代替 y=sin(t)./t; s=d*(cumtrapz(y)); %计算累计积分 plot(t,s),xlabel('x'),ylabel('y'),hold on y1=interp1(t,s,4.5) %插值法计算y(4.5) plot(4.5,y1,'*'),hold off 运行结果： y1 = 1.6542 4. 设0)0(,1)0(,1)(2)(3)(22===+-dt dy y t y dt t dy dt t y d ，用数值法和符号法求5.0)(=t t y 。 （提示：注意ode45和 dsolve 的用法。） 数值法： 编写M 函数文件DyDt.m function ydot=DyDt(t,y) ydot=[y(2);3* y(2)-2*y(1)-1]; %解算微分方程 tspan=[0.5,1]; y0=[1;0]; [tt,yy]=ode45(@DyDt,tspan,y0); figure(1) plot(tt,yy(:,1)) xlabel('t'),title('x(t)')

数值分析报告上机题(matlab版)(东南大学)

数值分析上机报告

第一章 一、题目 精确值为)1 1123(21+--N N 。 1) 编制按从大到小的顺序11 131121222-+ ??+-+-=N S N ，计算S N 的通用程序。 2) 编制按从小到大的顺序1 21 1)1(111222-+ ??+--+-= N N S N ，计算S N 的通用程序。 3) 按两种顺序分别计算64210,10,10S S S ,并指出有效位数。（编制程序时用单精度） 4) 通过本次上机题，你明白了什么？ 二、通用程序

三、求解结果 四、结果分析 可以得出，算法对误差的传播又一定的影响，在计算时选一种好的算法可以使结果更为精确。从以上的结果可以看到从大到小的顺序导致大数吃小数的现象，容易产生较大的误差，求和运算从小数到大数算所得到的结果才比较准确。

第二章 一、题目 （1）给定初值0x 及容许误差ε，编制牛顿法解方程f(x)=0的通用程序。 （2）给定方程03 )(3 =-=x x x f ,易知其有三个根3,0,3321=*=*- =*x x x a) 由牛顿方法的局部收敛性可知存在,0>δ当),(0δδ+-∈x 时，Newton 迭代序列收 敛于根x 2*。试确定尽可能大的δ。 b)试取若干初始值，观察当),1(),1,(),,(),,1(),1,(0+∞+-----∞∈δδδδx 时Newton 序列的收敛性以及收敛于哪一个根。 （3）通过本上机题，你明白了什么？ 二、通用程序

1.运行search.m 文件 结果为： The maximum delta is 0.774597 即得最大的δ为0.774597，Newton 迭代序列收敛于根* 2x =0的最大区间为 （-0.774597，0.774597）。 2.运行Newton.m 文件 在区间(,1),(1,),(,),(,1),(1,)δδδδ-∞----++∞上各输入若干个数，计算结果如下： 区间(,1)-∞-上取-1000，-100，-50，-30，-10，-8，-7，-5，-3，-1.5

数值分析五个题目的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数值计算-习题

1. 用函数roots求方程的根 roots([1 -1 -1]) 2. ，在n个节点(n不要太大，如取5~11)上用分段线性、三次方、样条插值方法，计算m个插值点(m可取50~100)的函数值。(注，n取10，m取100) x=linspace(0,2*pi,10); y=sin(x); xi=linspace(0,2*pi,100); y1=interp1(x,y,xi); y2=interp1(x,y,xi,'spline'); y3=interp1(x,y,xi,'cublic'); 3. 测得某地大气压强随高度变化的一组数据如表3-11 所示，试用插值法和拟合法估算高度为0，100，200，300，......，2000米时的大气压强值。 表3-11 某地大气压强随高度变化数据 高度/m 0 300 600 1000 1500 2000 压强/Pa x=[0 300 600 1000 1500 2000]; y=[ ]; xi=0:100:2000;

y1=interp1(x,y,xi,'spline'); p=polyfit(x,y,3); y2=polyval(p,xi); 4. 利用梯形法和辛普森法求定积分的值： 梯形法： x=linspace(-3,3,200); y=exp(-x.^2/2)/(2*pi); I1=trapz(x,y) 辛普森法： I2=quad('exp(-x.^2/2)/(2*pi)',-3,3) 或者： t='exp(-x.^2/2)/(2*pi)'; I2=quad(t,-3,3) 5. 分别用矩形法、梯形法、辛普森法和牛顿-科茨 4 种方法来近似计算定积分。 矩形法： x=linspace(0,1,100); y=x./(x.^2+4);

《数值分析与算法》例题matlab源程序

《数值分析与算法》喻文健，清华大学出版社，例题matlab源程序 %P33 example 2-2 a=1.0; b=1.5; while b-a>eps x=a+(b-a)/2 if f_ex22(a)*f_ex22(x)>0 a=x; else b=x; end end %P36 example 2-4 x=1.5 for i=1:10 x=(x+2)^(1/4) end %P42 example 2-7 x=1.5 for i=1:10 x=(3*x^4+2)/(4*x^3-1) end %P54 example2-9 x01=1; x02=2; x0=[x01;x02]; for i=1:8 x01=x0(1,1); x02=x0(2,1); f01=x01+2*x02-2; f02=x01^2+4*x02^2-4; J01=[1,2]; J02=[2*x01,8*x02]; s0=inv([J01;J02])*[f01;f02]x1=x0-s0 x0=x1; end P94 example 3-14 A=[5 -1 -1;-1 3 -1;-1 -1 5]; A(:,1)=A(:,1)./sqrt(A(1,1)); A(2,2)=sqrt(A(2,2)-A(2,1)^2); A(3,2)=(A(3,2)-A(3,1)^2)/A(2,2); A(3,3)=sqrt(A(3,3)-A(3,1)^2-A(3, 2)^2); %%% P111, 6 n=12; a=hilb(n); x=ones(n,1); b=a*x; b=b+10^(-7)*ones(n,1); L=chol(a,'lower'); %Cholesky y=ones(n,1); y(1)=b(1)/L(1,1); for i=2:n y(i)=(b(i)-sum(L(i,1:i-1)*y(1:i-1,1)))/L(i,i);%% L lower end M=L'; x(n)=y(n)/M(n,n); for i=n-1:-1:1 x(i)=(y(i)-sum(M(i,i+1:n)*y(i+1: n,1)))/M(i,i);%% L' upper end r=b-a*x; delt_x=x-ones(n,1); norm_r=norm(r,inf) % norm of r norm_delt_x=norm(delt_x,inf) % norm of delt_x

数值分析作业MATLAB程序题

function [x,y1,y2]=my(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y1(1)=y0; y2(1)=y0; for n=1:length(x)-1 y1(n+1)=y1(n)+h*feval(dyfun,x(n),y1(n)); k1=feval(dyfun,x(n),y2(n)); k2=feval(dyfun,x(n)+h/2,y2(n)+h/2*k1); k3=feval(dyfun,x(n)+h/2,y2(n)+h/2*k2); k4=feval(dyfun,x(n)+h,y2(n)+h*k3); y2(n+1)=y2(n)+h/6*(k1+2*k3+2*k2+k4); end x=x';y1=y1';y2=y2'; x3=0:0.02:1;y3=1/3*exp(-50*x)+x.^2; plot(x,y1,'r',x,y2,'g',x3,y3,'b'); legend('Euler','Runge-Kutta','真值解') >> dyfun=inline('-50*y+50*x.^2+2*x'); [x,y1,y2]=my(dyfun,[0,1],1/3,0.02);[x,y1,y2] H=0.02时

ans = 0 0.3333 0.3333 0.0200 0 0.1254 0.0400 0.0012 0.0485 0.0600 0.0032 0.0212 0.0800 0.0060 0.0130 0.1000 0.0096 0.0125 0.1200 0.0140 0.0153 0.1400 0.0192 0.0200 0.1600 0.0252 0.0257 0.1800 0.0320 0.0325 0.2000 0.0396 0.0400 0.2200 0.0480 0.0484 0.2400 0.0572 0.0576 0.2600 0.0672 0.0676 0.2800 0.0780 0.0784 0.3000 0.0896 0.0900 0.3200 0.1020 0.1024 0.3400 0.1152 0.1156 0.3600 0.1292 0.1296 0.3800 0.1440 0.1444 0.4000 0.1596 0.1600 0.4200 0.1760 0.1764 0.4400 0.1932 0.1936 0.4600 0.2112 0.2116 0.4800 0.2300 0.2304 0.5000 0.2496 0.2500 0.5200 0.2700 0.2704 0.5400 0.2912 0.2916 0.5600 0.3132 0.3136 0.5800 0.3360 0.3364 0.6000 0.3596 0.3600 0.6200 0.3840 0.3844 0.6400 0.4092 0.4096 0.6600 0.4352 0.4356 0.6800 0.4620 0.4624 0.7000 0.4896 0.4900 0.7200 0.5180 0.5184 0.7400 0.5472 0.5476 0.7600 0.5772 0.5776

MATLAB 数值分析 第二次 50页第二章计算实习题2

第二次 50页第二章计算实习题2 clear all clc x=-1:.1:1; y=1./(1+25*x.^2); plot(x,y) for n=10:10:20 x1=linspace(-1,1,n) y1=interp1(x,y,x1) hold on plot(x1,y1,'r') x2=linspace(-1,1,n); a=polyfit(x,y,3);%polyfit:给定n+1个点将可以唯一确定一个n 阶多项式。利用命令polyfit 可容易确定多项式的系数 y2=polyval(a,x2)%polyval:用命令polyval 计算多项式的值 hold on plot(x2,y2,'g') end 图像： -1-0.8-0.6-0.4-0.200.20.40.60.81 -0.20 0.2 0.4 0.6 0.8 1 1.2

对于给定的函数2 1()125f x x = +在区间[1,1]-上取值，试求三次曲线拟合并画出拟合曲线方程。 8次拟合曲线; 程序：clear all clc x=-1:0.01:1; f=1./(1+25*x.^2); plot(x,f,'*') i=0:10; x1=-1+0.2.*i; f1=1./(1+25*x1.^2); a=polyfit(x1,f1,8) f2=polyval(a,x1); hold on plot(x1,f2,'r') 曲线系数(从高到低)：20.4662 0.0000 -43.6006 -0.0000 30.8171 0.0000 -8.5318 -0.0000 0.8880 -1-0.8-0.6-0.4-0.200.20.40.60.81 3次拟合曲线; clear all clc x=-1:0.01:1; f=1./(1+25*x.^2); plot(x,f,'*') i=0:10;

(完整版)MATLAB数值分析实例

2016-2017 第一学期数值分析 上机实验报告 姓名： xxx 学号： 20162……. 学院：土木工程学院 导师：……….. 联系电话：………….. 指导老师：………..

目录 第一题 (1) 1.1题目要求 (1) 1.2程序编写 (1) 1.3计算结果及分析 (2) 第二题 (4) 2.1题目要求 (4) 2.2程序编写 (4) 2.3计算结果及分析 (6) 第三题 (7) 3.1题目要求 (7) 3.2程序编写 (7) 3.3计算结果及分析 (8) 第四题 (9) 4.1题目要求 (9) 4.2程序编写 (9) 4.3计算结果及分析 (10) 第五题 (11) 5.1题目要求 (11) 5.2程序编写 (12) 5.3计算结果及分析 (13) 第六题 (17) 6.1题目要求 (17) 6.2程序编写 (17) 6.3计算结果及分析 (18) 6.4程序改进 (18)

第一题 选做的是第（1）小问。 1.1题目要求 编出不动点迭代法求根的程序；把x3+4x2?10=0写成至少四种x=g(x)的形式，取初值x0=1.5,进行不动点迭代求根，并比较收敛性及收敛速度。1.2程序编写

1.3计算结果及分析 ①第一种迭代公式：x=x3+4x2+x?10；matlab计算结果如下：（以下为命令行窗口的内容） 2；matlab计算结果如下： ②第二种迭代公式：x=√(10?x3)/4 （以下为命令窗口内容） 2；matlab计算结果如下： ③第三种迭代公式：x=√10 （以下为命令窗口内容）

?；matlab计算结果如下： ④第四种迭代公式：x=10(x2+4x) （以下为命令窗口内容） 上述4种迭代公式，1、4两种由于在x真实值附近|g`(x)|>1,不满足迭代局部收敛条件，所以迭代序列不收敛。对于2、3两种式子，由于在x真实值附近|g`(x)|<=L<1,满足迭代局部收敛条件,所以迭代序列收敛。对于2、3两迭代公式，由于L3

数值分析matlab程序实例

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 =