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

2、设有某实验数据如下：

（1）在MATLAB中作图观察离散点的结构，用多项式拟合的方法拟合一个合适的多项式函数；

（2）在MATLAB中作出离散点和拟合曲线图.

解（1）：首先观察离散点的结构，matlab中的程序如下，

x=[-3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 ];

>> y=[-3.99 -3.3011 -2.4161 -1.4293 -0.4597 0.37758 1 1.3776 1.5403 1.5707 1.5839 1.6989 2.01]; >> plot(x,y,'r*')

图形如下：

离散点近似如抛物线，所以用二次多项式拟合，所以matlab程序如下：

x=[-3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 ];

>> y=[-3.99 -3.3011 -2.4161 -1.4293 -0.4597 0.37758 1 1.3776 1.5403 1.5707 1.5839 1.6989 2.01]; >> s=polyfit(x,y,2);

>> p=poly2str(s,'t')

p =

-0.22214 t^2 + 1 t + 0.74384

（2）做出离散点与拟合曲线的程序如下：

x=[-3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 ];

>> f=[-3.99 -3.3011 -2.4161 -1.4293 -0.4597 0.37758 1 1.3776 1.5403 1.5707 1.5839 1.6989 2.01]; >> p=polyfit(x,f,2); >> y=polyval(p,x); >> plot(x,f,'r+',x,y,'k'); >> xlabel('x'); >> ylabel('y'); >> title('拟合');

>> gtext(datestr(today))

得出的离散点与拟合曲线图像如下：

3、在MATLAB 中用复合Simpson 公式编程计算 2

21

()I x x dx -=

--⎰

（要求积分精度为510- ）

解：matlab 程序如下：

[I,step] = jfSimpson('-x-x^2',-1,2,2,1.0e-5) I =

-4.5000

step =

4

附复合Simpson 程序如下：

function [I,step] = jfSimpson(f,a,b,type,eps) %type = 1 辛普森公式 %type = 2 复合辛普森公式

if (type==2 && nargin==4)

eps=1.0e-4; %缺省精度为0.0001 end

I=0;

switch type case 1,

I=((b-a)/6)*(subs(sym(f),findsym(sym(f)),a)+... 4*subs(sym(f),findsym(sym(f)),(a+b)/2)+... subs(sym(f),findsym(sym(f)),b)); step=1;

case 2, n=2;

h=(b-a)/2; I1=0;

I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h; while abs(I2-I1)>eps n=n+1;

h=(b-a)/n; I1=I2; I2=0;

for i=0:n-1 x=a+h*i; x1=x+h;

I2=I2+(h/6)*(subs(sym(f),findsym(sym(f)),x)+... 4*subs(sym(f),findsym(sym(f)),(x+x1)/2)+... subs(sym(f),findsym(sym(f)),x1)); end end I=I2; step=n; end

4、在MATLAB 中用四阶Runge-Kutta 法编程求解微分方程初值问题

()2

3(03)01dy y x x dx y ⎧=-++≤≤⎪⎨⎪=⎩

， 并在MATLAB 中画图比较方程的解析解与R-K 解的结果。

解：第1步：首先先把经典RK4子程序编出来如下，用RK4.m 保存。