MATLAB语言基础与应用(第二版)第5章 习题答案

第5章习题与答案
5.1用矩阵三角分解方法解方程组
1231231
23214453186920
x x x x x x x x x +-=⎧⎪
-+=⎨⎪+-=⎩ 解答：
>>A=[2 1 -1;4 -1 3;6 9 -1] A =
2 1 -1 4 -1
3 6 9 -1 >>b=[1
4 18 20]; b =
14 18 20 >> [L, U, P]=lu(A) L =
1.0000 0 0 0.6667 1.0000 0 0.3333 0.2857 1.0000 U =
6.0000 9.0000 -1.0000 0 -
7.0000 3.6667 0 0 -1.7143 P =
0 0 1 0 1 0 1 0 0 >> y=backsub(L,P*b’) y =
20.0000 4.6667 6.0000 >> x=backsub(U,y) x =
6.5000 -2.5000 -3.5000 5.2 Cholesky 分解方法解方程组
123121
33235223
3127
x x x x x x x ++=⎧⎪
+=⎨⎪+=⎩ 解答：
>> A=[3 2 3;2 2 0;3 0 12] A =
3 2 3
2 2 0
3 0 12
>> b=[5;3;7]
b =
5
3
7
>> L=chol(A)
L =
1.7321 1.1547 1.7321
0 0.8165 -2.4495
0 0 1.7321
>> y=backsub(L,b)
y =
-11.6871 15.7986 4.0415
>> x=backsub(L',y)
x =
-6.7475 28.8917 49.9399
5.3
解答:
观察数据点图形
>> x=0:0.5:2.5
x =
0 0.5000 1.0000 1.5000 2.0000 2.5000 >> y=[2.0 1.1 0.9 0.6 0.4 0.3]
y =
2.0000 1.1000 0.9000 0.6000 0.4000 0.3000 >> plot(x,y)
图5.1 离
散点分布示意图
从图5.1观察数据点分布，用二次曲线拟合。

>> x=0:0.5:2.5 x =
0 0.5000 1.0000 1.5000 2.0000 2.5000 >> y=[2.0 1.1 0.9 0.6 0.4 0.3] y =
2.0000 1.1000 0.9000 0.6000 0.4000 0.3000 >> plot(x,y)
>> orthpolyfit(x,y,2) ans =
0.2857 -1.3371 1.9000
于是得到的二次拟合曲线方程为2
0.2857 1.3371 1.9y x x =-+ >> xx=-1:0.01:3;
>> pp=orthpolyfit(x,y,2); >> yy=polyval(pp,xx); >> plot(xx,yy) >> hold on
>> scatter(x,y,'k*')
00.51 1.52 2.5
图5.2 二次曲线拟合示意图
5.4已知函数()y f x =的观测数据为
试构造不超过三次的插值多项式，并计算(1)f -的近似值。

解答：
使用MATLAB 工具箱
-1
-0.500.51 1.52
2.53
00.511.522.533.54
图5.3 插值工具箱示意图
Linear model Poly3:
f(x) = p1*x^3 + p2*x^2 + p3*x + p4 Coefficients:
p1 = 0.119 p2 = -0.07143 p3 = -2.619 p4 = 1 >> p=[p1 p2 p3 p4]; >> polyval(p,-1) ans = 3.4286
即(1)f -=3.4286.
5.5用最小二乘方法对下表中的数据用经验公式bx
y ae =进行拟合
解答：
经验公式bx
y ae =首先线性化，令ln ,ln y y y a bx ==+则。

图5.4（a ）数据点分布图
对数据点{},i i x y 进行分析，
图5.4（b ）线性化数据点分布图
对图5.4（b ）数据点进行线性拟合 >> c=linearfit(x,yp,1) c =
12345678910
-0.3011 8.1694
于是ln 8.1694,0.3011a b ==-，求得3531.39,0.3011a b ==-。

于是所求的函数为0.30113531.39bx
x
y ae e
-==，拟合示意图如下。

图5.4（c ）数据点拟合示意图
5.6使用数值方法求函数22
1
()(1)f x x =+在 1.1,1.2,1.3x =处的导数值。

()f x 的函数值如
下表中：
在MA TLAB 命令行窗口中输入 >> x=[1.0 1.1 1.2 1.3 1.4];
>> y=[0.2500 0.2268 0.2066 0.1890 0.1736];
>>[A,df1]=npointdiff(x,y,1.1);[A,df2]=npointdiff(x,y,1.2);[A,df3]=npointdiff(x,y,1.3);[df1;df2;df3] ans =
-0.2320 -0.2170 -0.2033
5.7用Romberg 方法计算积分
2
1
sin()
x dx x
⎰
，要求结果保留5位小数。

解答：
编写被积函数M 文件myfun07.m function f=myfun07(x)
0246810
12
f=sin(x)./x;
参考EXAMP5008，修改romberg函数如下
function[s,n]=romberg(f,a,b,Eps)
if nargin<4
Eps=1e-6;
end
s=1;
s0=0;
k=2;
t=[];
t(1,1)=(b-a)*(feval(f,a)+feval(f,b))/2;
while abs(s-s0)>Eps
h=(b-a)/2^(k-1);
w=0;
if h~=0
for i=1:2^(k-1)-1
w=w+feval(f,a+i*h);
end
t(k,1)=h*(feval(f,a)/2+w+feval(f,b)/2);
for j=2:k
for i=1:(k-j+1)
t(i,j)=(4^(j-1)*t(i+1,j-1)-t(i,j-1))/(4^(j-1)-1);
end
end
s=t(1,k);
s0=t(1,k-1);
k=k+1;
n=k;
else
s=s0
n=-k
end
end
在命令行中调用romberg函数:
>> y=romberg('myfun07',1,2);vpa(y,5)
ans = 0.65933
5.8求矩阵
1.0 1.00.5
1.0 1.00.25
0.50.25 2.0
⎡⎤
⎢⎥
⎢⎥
⎢⎥
⎣⎦
的全部特征值及其对应的特征向量。

解答：
在MA TLAB命令行窗口中输入
>> A=[1.0 1.0 0.5;1.0 1.0 0.25;0.5 0.25 2.0]; >> [V D] = eig(A)
V =
0.7212 0.4443 0.5315 -0.6863 0.5621 0.4615 -0.0937 -0.6976 0.7103 D =
-0.0166 0 0 0 1.4801 0
0 0 2.5365
5.9使用Euler 方法解初值问题2(01)
(0)1
y y xy x y '⎧=--≤≤⎨=⎩，计算结果保留6位小数。

解答：
编写右端向量函数M 文件myfun11.m function f=myfun10(x,y) f=-y-x.*y.*y;
利用EXAMP5014，在命令行中调用euler 函数
[x y]=euler('myfun10',0,1,1,0.1)；format long ；y=vpa(y,6) y =
[ 1.0, 0.9, 0.8019, 0.708849, 0.62289, 0.545081, 0.475718, 0.414567, 0.36108, 0.314542, 0.274183]
5.10求解二阶微分方程的初值问题
,(01)(0)0,(0)1y y x x y y ''-=≤≤⎧⎨
'==⎩
解答：
令y z '=，原问题化为一阶方程组,(01)
(0)0,(0)1y z z x y x y z '=⎧⎪
'=+≤≤⎨⎪==⎩
编写右端向量函数M 文件myfun11.m function y=myfun11(x,w) y=zeros(2,1); %y=w(1) %z=w(2) y(1)=w(2); y(2)=x+w(1);
在命令行中调用ode45函数求原函数及导函数:
[x,y]=ode45(@dy,[0 1],[0 1]);x=x',z=y';y=z(1,:),yy=z(2,:) x =
Columns 1 through 15
0 0.0001 0.0001 0.0002 0.0002 0.0005 0.0007 0.0010 0.0012 0.0025 0.0037 0.0050 0.0062 0.0125 0.0188 Columns 16 through 30
0.0251 0.0313 0.0563 0.0813 0.1063 0.1313 0.1563 0.1813
0.2063 0.2313 0.2563 0.2813 0.3063 0.3313 0.3563
Columns 31 through 45
0.3813 0.4063 0.4313 0.4563 0.4813 0.5063 0.5313 0.5563 0.5813 0.6063 0.6313 0.6563 0.6813 0.7063 0.7313
Columns 46 through 57
0.7563 0.7813 0.8063 0.8313 0.8563 0.8813 0.9063 0.9313 0.9485 0.9657 0.9828 1.0000
y =
Columns 1 through 15
0 0.0001 0.0001 0.0002 0.0002 0.0005 0.0007 0.0010 0.0012 0.0025 0.0037 0.0050 0.0062 0.0125 0.0188
Columns 16 through 30
0.0251 0.0314 0.0564 0.0815 0.1067 0.1321 0.1576 0.1833 0.2093 0.2355 0.2620 0.2888 0.3160 0.3435 0.3715
Columns 31 through 45
0.4000 0.4289 0.4584 0.4884 0.5190 0.5502 0.5821 0.6146 0.6480 0.6820 0.7169 0.7526 0.7893 0.8268 0.8653
Columns 46 through 57
0.9048 0.9453 0.9869 1.0296 1.0735 1.1186 1.1649 1.2126
1.2460 1.2802 1.3149 1.3504
yy =
Columns 1 through 15
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0002 1.0004
Columns 16 through 30
1.0006 1.0010 1.0032 1.0066 1.0113 1.0173 1.0245 1.0330 1.0427 1.0538 1.0661 1.0797 1.0946 1.1108 1.1283
Columns 31 through 45
1.1472 1.1674 1.1890 1.2119 1.2362 1.2619 1.2890 1.3176 1.3476 1.3791 1.4120 1.4465 1.4825 1.5200 1.5591
Columns 46 through 57
1.5999 1.6422 1.6862 1.7319 1.7793 1.8284 1.8793 1.9320 1.9692
2.0073 2.0463 2.0862
作图
plot(x,y,'k-',x,yy,'k*');
title('Solution of ODEs with ode45')
legend('primitive function','derived function')
图形结果如下：
图5.5 二阶微分方程求解示意图
5.11取步长0.1h =，用经典公式求解初值问题
2,(01)(0)1y xy x y '=≤≤⎧⎨=⎩
解答：
参考EXAMP5016，编写右端函数M 文件myfun12.m
function f=myfun12(x,y)
f=2*x*y;
修改函数rungekutta 如下
function [xx,yy]=rungekutta(f,a,b,y0,h)
% f is the function entered as a string 'f'
% a,b are the left and right end points
% ya is the initial condition
% h is the step
x=a:h:b;
n=length(x);
y=zeros(1,n);
y(1)=y0;
for i=1:n-1
k1=feval(f,x(i),y(i));
k2=feval(f,x(i)+h/2,y(i)+h/2*k1);
k3=feval(f,x(i)+h/2,y(i)+h/2*k2);
k4=feval(f,x(i+1),y(i)+h*k3);
y(i+1)=(k1+2*k2+2*k3+k4)*h/6+y(i); end
xx=x;
yy=y;
ey=exp(x.*x)
plot(xx,yy,'kp',xx,ey,'k*');
xlabel('Independent Variable x'),ylabel('Dependent Variable y')
title('Solution of ODE with Runge-Kutta method')
legend('numal solution','exact solution')
在命令行中调用rungekutta函数:
>> [x y]=rungekutta('myfun12',0,1,1,0.1)
ey =
1.0000 1.0101 1.0408 1.0942 1.1735 1.2840 1.4333 1.6323 1.8965
2.2479 2.7183
x =
0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000
y =
1.0000 1.0101 1.0408 1.0942 1.1735 1.2840 1.4333 1.6323 1.8965
2.2479 2.7183
图形结果如下：
图5.6 经典公式求解示意图。