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

F=[3.2015,2.2560,879.5,1835.9,2968.8,4136.2,5237.9,6152.7,...

6725.3,6848.3,6403.5,6824.7,7328.5,7857.6];

F1=interp1(T,F,X) %用线性插值方法插值

F1=interp1(T,F,X,'nearest') %用最近点插值方法插值

F1=interp1(T,F,X,'spline') %用3次样条插值方法插值

F1=interp1(T,F,X,'cubic') %用3次多项式插值方法插值

例6.9设z=x2+y2，对z函数在[0,1]×[0,2]区域内进行插值。

x=0:0.1:1;y=0:0.2:2;

[X,Y]=meshgrid(x,y); %产生自变量网格坐标

Z=X.^2+Y.^2; %求对应的函数值

interp2(x,y,Z,0.5,0.5) %在(0.5,0.5)点插值

interp2(x,y,Z,[0.5 0.6],0.4) %在(0.5,0.4)点和(0.6，0.4)点插值

interp2(x,y,Z,[0.5 0.6],[0.4 0.5]) %在(0.5,0.4)点和(0.6，0.5)点插值

%下一命令在(0.5,0.4),(0.6,0.4),(0.5,0.5)和(0.6,0.5)各点插值

interp2(x,y,Z,[0.5 0.6]',[0.4 0.5])

例6.10某实验对一根长10米的钢轨进行热源的温度传播测试。用x表示测量点(米)，用h表示测量时间(秒)，用T表示测得各点的温度(℃)，测量结果如表6.2所示。

试用3次多项式插值求出在一分钟内每隔10秒、钢轨每隔0.5米处的温度。

x=0:2.5:10;

h=[0:30:60]';

T=[95,14,0,0,0;88,48,32,12,6;67,64,54,48,41];

xi=[0:0.5:10];

hi=[0:10:60]';

temps=interp2(x,h,T,xi,hi,'cubic');

mesh(xi,hi,temps);

例6.11 用一个3次多项式在区间[0,2π]内逼近函数x sin 。

X=linspace(0,2*pi,50); Y=sin(X);

P=polyfit(X,Y,3) %得到3次多项式的系数和误差 X=linspace(0,2*pi,20); Y=sin(X); Y1=polyval(P,X) plot(X,Y,':o',X,Y1,'-*')

例6.12 设

54322

()352756()353

f x x x x x x

g x x x =-+-++=+-

(1) 求f (x )+g (x )、f (x )-g (x )。 (2) 求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) conv(f,g) %求f(x)*g(x)

[Q,r]=deconv(f,g) %求f(x)/g(x)，商式送Q ，余式送r 。

例6.13 求有理分式的导数。

542109632

3585

()10567100

x x x x f x x x x x x +-+-=+++-- 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.14 已知多项式x 4+8x 3-10，分别取x =1.2和一个2×3矩阵为自变量计算该多项式的值。

A=[1,8,0,0,-10]; %4次多项式系数 x=1.2; %取自变量为一数值 y1=polyval(A,x)

x=[-1,1.2,-1.4;2,-1.8,1.6] %给出一个矩阵x

y2=polyval(A,x) % 分别计算矩阵x 中各元素为自变量的多项式之值

例6.15 仍以多项式x 4+8x 3-10为例，取一个2×2矩阵为自变量分别用polyval 和polyvalm

计算该多项式的值。

A=[1,8,0,0,-10]; %多项式系数 x=[-1,1.2;2,-1.8] %给出一个矩阵x y1=polyval(A,x) %计算代数多项式的值 y2=polyvalm(A,x) %计算矩阵多项式的值

例6.16 求多项式x 4+8x 3-10的根。

A=[1,8,0,0,-10]; x=roots(A)

例6.17 已知

52.7543)(235+--+=x x x x x f

(1) 计算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(x)

例6.18 设x 由[0,2π]间均匀分布的10个点组成，求x sin 的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)

例6.19 设

255122)(623+++++-+=x x x x x x f

用不同的方法求函数f (x )的数值导数，并在同一个坐标系中做出f '(x )的图像。

f=inline('sqrt(x.^3+2*x.^2-x+12)+(x+5).^(1/6)+5*x+2');

g=inline('(3*x.^2+4*x-1)./sqrt(x.^3+2*x.^2-x+12)/2+1/6./(x+5).^(5/6)+5'); x=-3:0.01:3;

p=polyfit(x,f(x),5); %用5次多项式p 拟合f(x) dp=polyder(p); %对拟合多项式p 求导数dp dpx=polyval(dp,x); %求dp 在假设点的函数值 dx=diff(f([x,3.01]))/0.01; %直接对f(x)求数值导数

gx=g(x); %求函数f 的导函数g 在假设点的导数

plot(x,dpx,x,dx,'.',x,gx,'-'); %作图

例6.20 用两种不同的方法求：

2

1

e d x I x -=?

先建立一个函数文件ex.m ：

function ex=ex(x) ex=exp(-x.^2);

然后在MATLAB 命令窗口，输入命令：

format long

I=quad('ex',0,1) %注意函数名应加字符引号 I=quadl('ex',0,1)

例6.21 用trapz 函数计算：

2

1

e

d x I x -=?

X=0:0.01:1; Y=exp(-X.^2); trapz(X,Y)

例6.22 计算二重定积分

2

1

2

/2

212

e sin()d d x

I x y x y ---=+?

?

(1) 建立一个函数文件fxy.m ：

function f=fxy(x,y) global ki;

ki=ki+1; %ki 用于统计被积函数的调用次数 f=exp(-x.^2/2).*sin(x.^2+y); (2) 调用dblquad 函数求解。

global ki;ki=0;

I=dblquad('fxy',-2,2,-1,1) ki

例6.23 给定数学函数

x (t )=12sin(2π×10t +π/4)+5cos(2π×40t )

取N =128，试对t 从0~1秒采样，用FFT 作快速傅立叶变换，绘制相应的振幅-频率图。

N=128; %采样点数 T=1; %采样时间终点

t=linspace(0,T,N); %给出N 个采样时间ti(i=1:N) x=12*sin(2*pi*10*t+pi/4)+5*cos(2*pi*40*t); % 求各采样点样本值x dt=t(2)-t(1); %采样周期 f=1/dt; %采样频率(Hz)

X=fft(x); %计算x 的快速傅立叶变换X F=X(1:N/2+1); %F(k)=X(k)(k=1:N/2+1) f=f*(0:N/2)/N; %使频率轴f 从零开始 plot(f,abs(F),'-*') %绘制振幅-频率图 xlabel('Frequency'); ylabel('|F(k)|')

例6.24 用直接解法求解下列线性方程组。

??????

?=--+=-+-=+-=+-+0

4662975135243214324

214321x x x x x x x x 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]'; x=A\b

例6.25 用LU 分解求解例6.24中的线性方程组。

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) [L,U ,P]=lu(A); x=U\(L\P*b)

例6.26 用QR 分解求解例6.24中的线性方程组。

A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; b=[13,-9,6,0]'; [Q,R]=qr(A); x=R\(Q\b) [Q,R,E]=qr(A); x=E*(R\(Q\b))

例6.27 用Cholesky 分解求解例6.24中的线性方程组。

A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; b=[13,-9,6,0]'; R=chol(A)

??? Error using ==> chol Matrix must be positive definite

Jacobi 迭代法的MATLAB 函数文件Jacobi.m 如下：

function [y,n]=jacobi(A,b,x0,eps) if nargin==3 eps=1.0e-6; elseif nargin<3 error return end

D=diag(diag(A)); %求A 的对角矩阵 L=-tril(A,-1); %求A 的下三角阵 U=-triu(A,1); %求A 的上三角阵 B=D\(L+U); f=D\b; y=B*x0+f;

n=1; %迭代次数 while norm(y-x0)>=eps x0=y; y=B*x0+f; n=n+1; end

例6.28 用Jacobi 迭代法求解下列线性方程组。设迭代初值为0，迭代精度为10-6。

??

?

??=+-=-+-=-6

10272109103232121x x x x x x x 在命令中调用函数文件Jacobi.m ，命令如下：

A=[10,-1,0;-1,10,-2;0,-2,10];

b=[9,7,6]';

[x,n]=jacobi(A,b,[0,0,0]',1.0e-6)

Gauss-Serdel 迭代法的MATLAB 函数文件gauseidel.m 如下：

function [y,n]=gauseidel(A,b,x0,eps) if nargin==3 eps=1.0e-6; elseif nargin<3 error return end

D=diag(diag(A)); %求A 的对角矩阵 L=-tril(A,-1); %求A 的下三角阵 U=-triu(A,1); %求A 的上三角阵 G=(D-L)\U; f=(D-L)\b; y=G*x0+f;

n=1; %迭代次数 while norm(y-x0)>=eps x0=y; y=G*x0+f; n=n+1; end

例6.29 用Gauss-Serdel 迭代法求解下列线性方程组。设迭代初值为0，迭代精度为10-6。

??

?

??=+-=-+-=-6

10272109103232121x x x x x x x 在命令中调用函数文件gauseidel.m ，命令如下：

A=[10,-1,0;-1,10,-2;0,-2,10]; b=[9,7,6]';

[x,n]=gauseidel(A,b,[0,0,0]',1.0e-6)

例6.30 分别用Jacobi 迭代和Gauss-Serdel 迭代法求解下列线性方程组，看是否收敛。

????

?

?????=????????????????????-679122111221321x x x a=[1,2,-2;1,1,1;2,2,1]; b=[9;7;6];

[x,n]=jacobi(a,b,[0;0;0]) [x,n]=gauseidel(a,b,[0;0;0])

有了上面这些讨论，下面设计一个求解线性方程组的函数文件line_solution.m 。

function [x,y]=line_solution(A,b) [m,n]=size(A); y=[];

if norm(b)>0 %非齐次方程组

if rank(A)==rank([A,b])

if rank(A)==n %有惟一解

disp('原方程组有惟一解x'); x=A\b;

else %方程组有无穷多个解，基础解系

disp('原方程组有无穷个解，特解为x ，其齐次方程组的基础解

系为y');

x=A\b; y=null(A,'r'); end else

disp('方程组无解'); %方程组无解 x=[]; end

else %齐次方程组

disp('原方程组有零解x'); x=zeros(n,1); %0解 if rank(A)

disp('方程组有无穷个解，基础解系为y'); %非0解

y=null(A,'r'); end end

例6.31 求解方程组。

???

??=-++=-+-=-+-3

22223531

324321

43214321x x x x x x x x x x x x A=[1,-2,3,-1;3,-1,5,-3;2,1,2,-2]; b=[1;2;3];

[x,y]=line_solution(A,b)

例6.32 求方程组的通解。

???

??=--+=+--=--+0

89544331

34321

43214321x x x x x x x x x x x x format rat %指定有理式格式输出 A=[1,1,-3,-1;3,-1,-3,4;1,5,-9,-8]; b=[1,4,0]';

[x,y]=line_solution(A,b); x,y

format short %恢复默认的短格式输出

例6.33 求1

()5f x =x +x

-在x 0=-5和x 0=1作为迭代初值时的零点。

先建立函数文件fz.m ：

function f=fz(x) f=x-1/x+5;

然后调用fzero 函数求根。：

fzero('fz',-5) %以-5作为迭代初值 fzero('fz',1) %以1作为迭代初值

例6.34 求下列方程组在(1,1,1)附近的解并对结果进行验证。

2sin e 0

0x x y z x y z xyz ?++=?

++=??=?

首先建立函数文件myfun.m 。

function F=myfun(X) x=X(1); y=X(2); z=X(3);

F(1)=sin(x)+y+z^2*exp(x); F(2)=x+y+z; F(3)=x*y*z;

在给定的初值x 0=1，y 0=1，z 0=1下，调用fsolve 函数求方程的根。

X=fsolve('myfun',[1,1,1],optimset('Display', 'off'))

例6.35 求圆和直线的两个交点。

圆：9222=++z y x

直线：???=---=++0

1630

653z y x z y x

先建立方程组函数文件fxyz.m ：

function F=fxyz(X) x=X(1); y=X(2); z=X(3);

F(1)=x^2+y^2+z^2-9; F(2)=3*x+5*y+6*z; F(3)=x-3*y-6*z-1;

再在MATLAB 命令窗口，输入命令：

X1=fsolve('fxyz',[-1,1,-1],optimset('Display', 'off')) %求第一个交点 X2=fsolve('fxyz',[1,-1,1],optimset('Display', 'off')) %求第二个交点

例6.36 求函数

1()5f x x x

=-

+ 在区间(-10,1)和(1,10)上的最小值点。

首先建立函数文件fx.m ：

function f=f(x) f=x-1/x+5;

上述函数文件也可用一个语句函数代替：

f=inline('x-1/x+5')

再在MATLAB 命令窗口，输入命令：

fminbnd('fx',-10,-1) %求函数在(-10，-1)内的最小值点和最小值 fminbnd(f,1,10) %求函数在(1，10)内的最小值点。注意函数名f 不用

加' 例6.37 设

2

2

2

4(,,)y z x y z f x y z x =+++

求函数f 在(0.5,0.5,0.5)附近的最小值。

建立函数文件fxyz.m ：

function f=fxyz(u) x=u(1);y=u(2);z=u(3); f=x+y.^2./x/4+z.^2./y+2./z; 在MALAB 命令窗口，输入命令：

[U,fmin]=fminsearch('fxyz',[0.5,0.5,0.5]) %求函数的最小值点和最小值

例6.38 求解有约束最优化问题。

3

1212

22120

,05

.05.04

.05.0..30

14.0)(min

21

2121x x x x x x x f x x x x x x t s x +

-++=???

??≥≥≥+≥+ 首先编写目标函数M 文件fop.m 。

function f=fop(x)

f=0.4*x(2)+x(1)^2+x(2)^2-x(1)*x(2)+1/30*x(1)^3;

再设定约束条件，并调用fmincon 函数求解此约束最优化问题。

x0=[0.5;0.5]; A=[-1,-0.5;-0.5,-1]; b=[-0.4;-0.5]; lb=[0;0];

option=optimset; https://www.wendangku.net/doc/1d10136128.html,rgeScale='off'; option.Display ='off'; [x,f]=fmincon('fop ',x0,A,b,[],[],lb,[],[],option)

例6.39 设有初值问题：

22' 014(1)y t y t t --=≤≤+，

y (0)=2

试求其数值解，并与精确解相比较(精确解为y (t )=11++t )。

(1) 建立函数文件funt.m 。

function y=funt(t,y) y=(y^2-t-2)/4/(t+1); (2) 求解微分方程。

t0=0;tf=10; y0=2;

[t,y]=ode23('funt',[t0,tf],y0); %求数值解 y1=sqrt(t+1)+1; %求精确解 plot(t,y,'b.',t,y1,'r-'); %通过图形来比较

例6.40 已知一个二阶线性系统的微分方程为：

22d 0,0

d (0)0'(0)1x

ax a t

x x ?+=>???==?， 其中a =2，绘制系统的时间响应曲线和相平面图。

建立一个函数文件sys.m ：

function xdot=sys(t,x) xdot=[-2*x(2);x(1)]; 取t0=0，tf=20，求微分方程的解：

t0=0;tf=20;

[t,x]=ode45('sys',[t0,tf],[1,0]); [t,x]

subplot(1,2,1); plot(t,x(:,2)); %解的曲线，即t-x subplot(1,2,2); plot(x(:,2),x(:,1)) %相平面曲线，即x-x’ axis equal

例6.41 设

200000000

00

0050010010

000

5????????=??-????-??

X 将X 转化为稀疏存储方式。

X=[2,0,0,0,0;0,0,0,0,0;0,0,0,5,0;0,1,0,0,-1;0,0,0,0,-5];

A=sparse(X)

例6.42 根据表示稀疏矩阵的矩阵A ，产生一个稀疏存储方式矩阵B 。

2213114335386612????-????=????????

A

A=[2,2,1;3,1,-1;4,3,3;5,3,8;6,6,12]; B=spconvert(A)

例6.43 求下列三对角线性代数方程组的解。

??

?

??????

???????=????????????????????????????????51230112624611413254321x x x x x B=[1,2,0;1,4,3;2,6,1;1,6,4;0,1,2]; %产生非0对角元素矩阵 d=[-1;0;1]; %产生非0对角元素位置向量 A=spdiags(B,d,5,5) %产生稀疏存储的系数矩阵 f=[0;3;2;1;5]; %方程右边参数向量 x=(inv(A)*f)' %求解

数值计算方法实验指导(Matlab版)

《数值计算方法》实验指导 （Matlab 版） 肇庆学院数学与统计学学院 计算方法课程组

1. 实验名称 实验1 算法设计原则验证(之相近数相减、大数吃小数和简化计算步骤) 2. 实验题目 有效数字的损失. 123 )与1000个较小的数(3 10 15)的和，验证 大数吃小数的现象. (3)分别用直接法和秦九韶算法计算多项式 P(x) a 0x n a 1x n 1 在x =1.00037 处的值?验证简化计算步骤能减少运算时间. n 1 对于第(3)题中的多项式P (x )，直接逐项计算需要n (n 1) 2 1 次乘法 和n 次加法，使用秦九韶算法 P(x) (((a °x ajx a 2)x a . 则只需要n 次乘法和n 次加法. 3. 实验目的 验证数值算法需遵循的若干规则. 4. 基础理论 设计数值算法时，应避免两个相近的数相减、防止大数吃小数、简化计算步骤减少运算 次数以减少运算时间并降低舍入误差的积累. 两相近的数相减会损失有效数字的个数， 用一 《数值计算方法》实验 1报告 班级： 20xx 级 XXXXx 班 学号： 20xx2409xxxx 姓名： XXX 成绩: ⑴取 z 1016，计算z 1 Z 和 1/（、z 1 Z)，验证两个相近的数相减会造成 (2)按不同顺序求一个较大的数( a n 1 X a n

个大数依次加小数，小数会被大数吃掉，乘法运算次数太多会增加运算时间. 5.实验环境 操作系统：Win dows xp ；程序设计语言：Matlab 6.实验过程 (1)直接计算并比较； (2)法1 :大数逐个加1000个小数，法2 :先把1000个小数相加再与大数加； (3)将由高次项到低次项的系数保存到数组A[n]中，其中n为多项式次数. 7.结果与分析 (1)计算的~1V Z = _______________________________ ,1/( ~1 < z) ____________________ . 分析： (2)123逐次加1000个3 10 6的和是_________________________ ，先将1000个3 10 6相 加，再用这个和与123相加得_______________________ . 分析： (3)计算__________ 次的多项式： 直接计算的结果是___________________ ，用时___________________ ; 用秦九韶算法计算的结果是____________________ ，用时 ___________________ 分析：

数值分析Matlab作业

数值分析编程作业

2012年12月 第二章 14.考虑梯形电阻电路的设计，电路如下： 电路中的各个电流{i1,i2,…,i8}须满足下列线性方程组： 12 123 234 345 456 567 678 78 22/ 2520 2520 2520 2520 2520 2520 250 i i V R i i i i i i i i i i i i i i i i i i i i -= -+-= -+-= -+-= -+-= -+-= -+-= -+= 这是一个三对角方程组。设V=220V，R=27Ω，运用追赶法，求各段电路的电流量。Matlab程序如下： function chase () %追赶法求梯形电路中各段的电流量 a=input('请输入下主对角线向量a='); b=input('请输入主对角线向量b='); c=input('请输入上主对角线向量c='); d=input('请输入右端向量d='); n=input('请输入系数矩阵维数n='); u(1)=b(1); for i=2:n l(i)=a(i)/u(i-1); u(i)=b(i)-c(i-1)*l(i); end y(1)=d(1); for i=2:n y(i)=d(i)-l(i)*y(i-1); end x(n)=y(n)/u(n); i=n-1; while i>0 x(i)=(y(i)-c(i)*x(i+1))/u(i); i=i-1; end x 输入如下:

请输入下主对角线向量a=[0,-2,-2,-2,-2,-2,-2,-2]; 请输入主对角线向量b=[2,5,5,5,5,5,5,5]; 请输入上主对角线向量c=[-2,-2,-2,-2,-2,-2,-2,0]; 请输入方程组右端向量d=[220/27,0,0,0,0,0,0,0]; 请输入系数矩阵阶数n=8 运行结果如下： x = 8.1478 4.0737 2.0365 1.0175 0.5073 0.2506 0.1194 0.0477 第三章 14.试分别用(1)Jacobi 迭代法；(2)Gauss-Seidel 迭代法解线性方程组 1234510123412191232721735143231211743511512x x x x x ?????? ??????---????????????=--?????? --?????? ??????---?????? 迭代初始向量 (0)(0,0,0,0,0)T x =。 （1）雅可比迭代法程序如下： function jacobi() %Jacobi 迭代法 a=input('请输入系数矩阵a='); b=input('请输入右端向量b='); x0=input('请输入初始向量x0='); n=input('请输入系数矩阵阶数n='); er=input('请输入允许误差er='); N=input('请输入最大迭代次数N='); for i=1:n for j=1:n if i==j d(i,j)=a(i,j); else d(i,j)=0; end end end m=eye(5)-d\a; %迭代矩阵 g=d\b; x=m*x0+g; k=1; while k<=N %进行迭代 for i=1:5 if max(abs(x(i)-x0(i))) >er x=m*x+g; k=k+1;

数值分析的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 *

第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程序

列主元法 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 编程题库 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样卷答案

腹有诗书气自华 《数值计算方法与Matlab 》 样卷答案 一．填空题：（每空3分,共42分） 1． 8，6105.0-? 。 2．)(3)1(2)1(1)(3)1(2)1(1)(3)1(3)(3)(2)1(1)(3)(2)1(1)(2)1(2)(3)(2)(1)(3)(2)(1)(1)1(1)1(22)22()1()1(222)1()222(k k k k k k k k k k k k k k k k k k k k k k k k x x x x x x x x x x x x x x x x x x x x x x x x ωωωωωωωωωω ωωωω-+--=---?+=+--+-=---?+=++--=+--?+=+++++++++， )2,1(∈ω。 3．],[1b a C S m -∈。4. 1e 2e ---x ，???==-=?--? ,3,2,1,0;0,e 1d )(e 110k k x x g k x ，正交投影。 5. 2阶，6阶。 6．10.6658，10.9521，10.9501。 7. 4002.2)00.1(=ε，4030.2)01.1(=ε。 二．解下列各题：（每题9分，共36分） 1．解：令)1(2 3+=t x ， （2分） 则??-+++=+1123 02 dt )1(25.21)1(49d 1t t x x x ???++++???++-+-≈22)6.01(25.21)6.01(9525.219 8)6.01(25.21)6.01(9549 （8分） 210631.10≈ （9分） 2．解：记系数矩阵为A, 对增广矩阵[]b A |作初等行运算, ??????????--401533933112??????????--==5.55.115 .35.405.75.401125.1,5.11,31,2l l ??????????---=45.114005.75.4011212,3l , 所以13-=x ,2)5.75.1(5.4112=-=x x ,1)1(2 1321=-+-=x x x ,即方程组的解为 [1,2,-1]T . （4分） 故系数矩阵A 的LU 分解为???? ??????--???????????---=4005.75.40112115.1015.1001A 。 （6分）

数值分析 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/1d10136128.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 版)

《数值计算方法》 实验指导 （Matlab版） 肇庆学院数学与统计学学院 计算方法课程组

《数值计算方法》实验1报告 班级： 20xx 级XXXXx 班 学 号： 20xx2409xxxx 姓名： XXX 成绩： 1. 实验名称 实验1 算法设计原则验证（之相近数相减、大数吃小数和简化计算步骤） 2. 实验题目 (1) 取16 10=z ，计算 z z -+1和) 1/( 1z z ++，验证 两个相近的数相减会造成有效数字的损失． (2) 按不同顺序求一个较大的数（123）与1000个较小的数（15 310-?）的和，验证大数吃小 数的现象． (3) 分别用直接法和秦九韶算法计算多项式 n n n n a x a x a x a x P ++++=--1 1 1 )(Λ 在x =1.00037处的值．验证简化计算步骤能减少运算时间． 对于第(3)题中的多项式P (x )，直接逐项计算 需要21 12)1(+=+++-+n n n Λ次乘法和n 次加法，使用秦

九韶算法 n n a x a x a x a x a x P ++++=-)))((()(1210Λ 则只需要n 次乘法和n 次加法． 3. 实验目的 验证数值算法需遵循的若干规则． 4. 基础理论 设计数值算法时，应避免两个相近的数相减、防止大数吃小数、简化计算步骤减少运算次数以减少运算时间并降低舍入误差的积累．两相近的数相减会损失有效数字的个数，用一个大数依次加小数，小数会被大数吃掉，乘法运算次数太多会增加运算时间． 5. 实验环境 操作系统：Windows xp ； 程序设计语言：Matlab 6. 实验过程 (1) 直接计算并比较； (2) 法1：大数逐个加1000个小数，法2：先把1000个小数相加再与大数加； (3) 将由高次项到低次项的系数保存到数组A[n]中，其中n 为多项式次数．

解非线性方程的数值计算方法 用Matlab实现

湖北民族学院《数值分析》实验报告 实验名称 解非线性方程的数值计算方法 实验时间 2011年 11月 9日 姓名 王亚敏 班级 0209408 学号 020940807 成绩 实验报告内容要求： 一、实验目的与要求；二、实验内容；三、算法描述(数学原理或设计思路、计算公式、计算步骤)； 四、程序代码；五、数值结果；六、计算结果分析（如初值对结果的影响；不同方法的比较；该方法的特点和改进等）；七、实验中出现的问题，解决方法及体会（整个实验过程中（包括程序编写，上机调试等）出现的问题及其处理等广泛的问题）. 一、实验目的： ① 掌握二分法、不动点迭代法、牛顿迭代法、Steffensen 迭代法等常用的非线性方程迭代算法； ② 培养编程与上机调试能力. 二、实验内容： 1．请分别用二分法、不动点迭代法、牛顿迭代法、Steffensen 迭代法求方程 32 ()330f x x x x =+--=在1.5 附近的根. 参考答案：原方程在1.5 附近的根为1.732051x = 2．请自选两个非线性方程，分别用二分法、不动点迭代法、牛顿迭代法、Steffensen 迭代法求解. 三、算法 （一） 1．二分法 计算()0f x =的二分法如下： ① 输入求根取间[,]a b 和误差控制量ε，定义函数. 如果 ()()0f a f b <，转②；否则退出选用其它求根方法 ② 当a b ε->时，计算中点()/2x a b =+以及的值； 分情况处理 ()f x ε<：停止计算，，转④ ()()0f a f x <：修正区间[,][,]ax ab → ()()0f x f b <：修正区间[,][,]xb ab → ③ *2a b x += ④ 输出近似根* x 2．不动点(简单)迭代法算法思想 先将f(x)=0化成x=φ(x),然后取x 0,计算 x n+1=φ(x n ) n=0,1,2,… 若要求x * 满足f(x * )=0,则x * =)(*x ?，称x * 为函数)(x ?的一个不动点。

数值计算方法matlab程序

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> fun1=inline('(x+1)^(1/3)'); >> [x0,k]=iterate1(fun1,1.5) x0 = 1.3247 k = 7 >> fun1=inline('x^3-1'); >> [x0,k]=iterate1(fun1,1.5) x0 = Inf k = 9 >> Steffesen加速迭代（简单迭代法的加速）function [x0,k]=steffesen1(fun1,x0,ep,N) if nargin<4 N=500; end if nargin<3 ep=1e-5; end x=x0; x0=x+2*ep; k=0; while abs(x-x0)>ep & k

数值分析重要算法的matlab程序

数值分析重要算法的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)); 追赶法解三对角方程组

数值计算方法matlab 第二章 求根

1 第二章作业 问题描述： 不同温度的两种流体进入混合器混合，流出时具有相同的温度。流体A 和B 的热容（单位：cal/(mol ·K)）分别为： 2623.381 1.80410 4.30010pA c T T --=+?-? 1528.592 1.29010 4.07810pB c T T --=+?-? 焓变（单位：cal/mol ）为2 1 T p T H c dT ?= ? 。 A 进入混合器的温度为400℃， B 进入混合器的温度为700℃，A 的量（mol ）是B 的量（mol ）的两倍，试确定流体离开混合器的温度。 问题分析： 初始情况下，气体A 的温度比气体B 的温度低，故在混合过程中，气体A 温度升高，气体B 温度降低。由于没有外界加热或者做功，混合气体整体的焓变应该为零。 设A 有2mol ，B 有1mol ，根据焓变公式计算得到： 2 1 -262400 -22632= 6.762+3.608108.60010)6.762 1.80410 2.867105407.712T T A pA T H c dT T T dT T T T --?=?-?=+?-?-??（ 2 1 -152700 -1253=+1.29010 4.07810)0.64510 1.3591032958.030 T T B pB T H c dT T T dT T T T --?=?-?=+?-?-??（8.5928.592 而0A B H H ?+?=，故该问题最后变成求解方程 2263()15.3548.2541016.4571038365.742f T T T T --=+?-?- 的根的问题。接下来将采用二分法、试位法以及牛顿法进行改方程的求解。方程保存为f.m ，可在压缩文件中找到。 一、 二分法 二分法的基本思想为，确定有根区间，然后不断将区间二等分，通过判断f(x)的符号，逐步将区间缩小，直到有根区间足够小，便可满足精度要求的近似根。 本例中，可以清楚的得到有根区间为(400,700)。取容限误差为-3 0.510%?，可以保证5 位有效数字。程序编写存储于bisec.m 。 其中，bisec 函数定义为： function bisec(f_name,a,c,xmin,xmax,n_points) 调用时： >> bisec('f',400,700,400,700,1000) 相当于取了a=400;c=700;作图时横坐标取得是从400~700的范围，采样点为1000个。

数值分析上机题(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

数值计算方法与Matlab样卷答案

《数值计算方法与Matlab 》 样卷答案 一．填空题：（每空3分,共42分） 1． 8，6105.0-? 。 2．) (3 )1(2)1(1)(3)1(2)1(1)(3)1(3) (3)(2)1(1)(3)(2)1(1)(2) 1(2 ) (3)(2)(1)(3)(2)(1)(1)1(1)1(22)22()1()1(222)1()222(k k k k k k k k k k k k k k k k k k k k k k k k x x x x x x x x x x x x x x x x x x x x x x x x ωωωωω ωωωωω ωωωω-+--=---?+=+--+-=---?+=++--=+--?+=+++++++++， )2,1(∈ω。 3．],[1 b a C S m -∈。4. 1e 2e --- x ，???==-=?--? ,3,2,1,0;0,e 1d )(e 110 k k x x g k x ，正交投影。 5. 2阶，6阶。 6．10.6658，10.9521，10.9501。 7. 4002.2)00.1(=ε，4030.2)01.1(=ε。 二．解下列各题：（每题9分，共36分） 1．解：令)1(2 3 += t x ， （2分） 则 ?? -+++=+1 1 230 2 dt )1(25.21)1(49d 1t t x x x ? ??++++???++ -+-≈ 22 )6.01(25.21)6.01(9525.2198)6.01(25.21)6.01(9549 （8分） 210631.10≈ （9分） 2．解：记系数矩阵为A, 对增广矩阵[]b A |作初等行运算, ???????? ??--4015 33933112??????? ???--==5.55.115 .35.405 .75.401 125.1,5.11,31,2l l ???? ??????---=45.114005.75.401 1212,3l , 所以13-=x ,2)5.75.1(5.41 12=-=x x ,1)1(2 1 321=-+-=x x x ,即方程组的解为 [1,2,-1]T . （4分） 故系数矩阵A 的LU 分解为???? ??????--???????????---=4005.75.40112115.1015.1001 A 。 （6分） 由于∞ ∞ ∞-∞∞∞∞∞∞??=?≤?||||||||||||||||||||||||)(cond ||||||||1b b A A b b A x x ，

数值计算方法 Matlab实题训练(内附程序,模型)

《数值计算方法训练》 实习报告 题目：6－A组 院系：上海电力学院数理学院 专业年级：信息与计算科学专业2009级 学生姓名：XX远学号：20092426 2011年7月8日

第1题：含炭量与时间的关系 在某冶炼过程中，钢的含炭量y 与时间t 的统计数据如下 t 0 5 10 15 20 25 30 35 40 45 50 55 y 0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.5 8 4.02 4.64 （1）画出原始数据分布趋势图； （2）用最小二乘法求钢的含炭量y 与时间t 的拟合曲线32ct bt at y ++=； （3）打印出拟合曲线； （4）另外选用b at y =进行拟合，比较二种拟合的效果。 解：分析：使用到曲线拟合的最小二乘法，对于拟合函数，尽量转化为可以方便提炼出基函数的方程。在明确基函数的基础上，通过计算，得到各个系数，得到法方程组 （1），程序：function yuan(y) t=[0:5:55]； plot(t,y,'*') legend('原始数据分布趋势图') 运行结果：yuan([0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64]) 图1 原始数据分布趋势 （2），使用最小二乘法，就必须先取基函数，对于该题流程如下：

①：取基函数为：t =0? 21t =? 32t =? ②：由基函数和y 求法方程组的系数： ?????? ??? ? ??? ?=?=?=?=?=?=?=?=?=∑∑∑∑∑∑∑∑∑=========)(),()(),()(),() ()(),()()(),()()(),() ()(),()()(),()()(),(212121121101210212 1222212 11211121111212 1 02011210100121000i i i i i i i i i i i i i i i i i i i i i i i i i i i x y f x y f x y f x x x x x x x x x x x x ?????????????????????????????? ③：由这些系数，确定法方程组： B X A =??? ?? ? ????? ???? ???????=???? ? ?????=),(),(),(),(),(),(),(),(),(),(),(),(210222120121110020100?????????????????????f f f B A ④：解这个法方程组： ?? ?? ? ?????==?c b a X B X A ，得到拟合函数：3 2ct bt at y ++= 程序：function [a,b,c]=xian(y0) t0=[0:5:55]; k1=t0; k2=t0.*t0; k3=t0.*t0.*t0; A=[sum(k1.*k1) sum(k2.*k1) sum(k3.*k1);sum(k1.*k2) sum(k2.*k2) sum(k3.*k2);sum(k1.*k3) sum(k2.*k3) sum(k3.*k3)]; B=[sum(k1.*y0);sum(k2.*y0);sum(k3.*y0)]; x=pinv(A)*B; a=x(1,1); b=x(2,1); c=x(3,1); t=0:55; y=a.*t+b.*t.^2+c.*t.^3; plot(t,y,'--') hold on plot(t0,y0,'*')

《数值计算方法》(第二版) 数值计算实验 matlab 有截图

数值计算方法 （第二版） 数值实验 学院信息学院专业计算机应用技术姓名张书涵学号12077503005 任课老师孙云平使用软件MA TLAB 7.0 2012 年12 月15日

一、实验目的和意义 1、通过上机学会matlab数学软件的使用，根据算法2.2的思想和步骤，实现解三对角方程组的追赶法。 2、运用所学的三种迭代法即Jacobi迭代法、Gauss-Seidel迭代法、共轭梯度法算法的思想和步骤，解各类线性方程组，编出算法程序。 3、根据插值法的思想，掌握样条插值的步骤和思想，并根据样本值用matlab软件绘制图形，熟悉了matlab的绘图函数。 4、只有通过动手编程，才能意识到自己对课本上知识的欠缺。通过本次实验，复习本学期所学知识，查漏补缺。 二、实验环境 matlab7.0 三、实验结果 第一章误差 1.利用循环语句，计算数列√5，√√5，√√√5，……的极限，要求误差小于10^-8。程序设计如下： x =sqrt(5); m=0; while abs(x-m)>1e-8 m=result; x=sqrt(x); end x 运行结果：

第二章 解线性方程组的直接解法 14. 考虑梯形电阻电路的设计，电路如下： 电路中的各个电流{1i ，2i ，…，8i }须满足下列线性方程组： R V i i =- 22 21 0 252321=-+-i i i 0 252 432=-+-i i i 0 252 543=-+-i i i 0 252 654=-+-i i i 0 252 765=-+-i i i 0 252 876=-+-i i i 052 87=+-i i 这是一个三对角方程组。设V 220=V ，Ω=27R ，运用追赶法，求各段电路的电流量。 程序设计如下： a=[0 -2 -2 -2 -2 -2 -2 -2]; b=[2 5 5 5 5 5 5 5]; c=[-2 -2 -2 -2 -2 -2 -2]; d=[8.1481 0 0 0 0 0 0 0]; for i=2:8 a(i)=a(i)/b(i-1); b(i)=b(i)-c(i-1)*a(i); d(i)=d(i)-a(i)*d(i-1); end; d(8)=d(8)/b(8); for i=7:-1:1 d(i)=(d(i)-c(i)*d(i+1))/b(i); end;