>> [x,y]=meshgrid(-2:0.1:2,-2:0.1:2); %产生网格点 >> z=x.*exp(-x.^2-y.^2); >> plot3(x,y,z); 例 3:
>> x=-8:0.5:8; y=x'; >> a=ones(size(y))*x; >> b=y*ones(size(x)); >> c=sqrt(a.^2+b.^2)+eps; >> z=sin(c)./c; >> surf(x,y,z) >> mesh(x,y,z) 3) 特殊的三维图形函数
* [x,y,z]=sphere(n) % 画球,n默认值20 例1 : >> [a,b,c]=sphere(40); >> surf(a,b,c) >> axis('equal'); >> axis('square');
* [x,y,z]=cylinder(R,N) %R 母线向量,N分格数 例2 : >> x=0 :pi/20 :pi*3; >> r=5+cos(x);
>> [a,b,c]=cylinder(r,30); >> mesh(a,b,c) 练习：
例：需要考察函数sin（x）与它的各阶Taylor展开式
x3x5 sin(x)?x???...的近似情况。请编写一个函数（程序）。当输入Taylor展
3!5!x3
开式的阶数后，函数将自动画出y=sin(x)，y=x, y=x?，直到要求的次数位置的
3!
各阶展开式在[??,?]上的图形.
程序代码:
function plot_sin_taylor(n) x=-pi:pi/20:pi; plot(x,sin(x));
6
hold on
m=fix((n+1)/2); y=0; for k=1:m
y=y+(-1)^(k-1)/prod([1:2*k-1])*x.^(2*k-1); plot(x,y); end
axis([-pi pi -1.5 1.5]); xlabel('x');ylabel('y'); hold off
作业:
xe请仿照上例考察与它的各阶Taylor展开式在[-1,1]内的近似情况.
2nxxxe?1?x??...?
2!n!
数值分析实验二（上）
第三章 解线性方程组的直接方法
1. 实验目的
1)熟悉用Matlab 向量运算 2)用Matlab矩阵求逆及三角分解 3)掌握解线性方程组的直接方法
2. 试验内容
(1) Gauss_Jordan列主元消去法求方阵逆.
选择列主元, 消去对角线上方和下方的元素. 用其实现矩阵求逆如下: function [A,det]=GaJo_inv(A)
% Gauss-Jordan列主元消去法求方阵逆 P178 % [A,det]=GaJo_inv(A)
% A 要求逆矩阵
7
% det 按需求返回 A 的行列式 % A 返回逆矩阵,放在A 中 n=size(A); if n(1)~=n(2) error('不是方阵!'); end
n=n(1);det=1; flag=1:n; for k=1:n
t=find(abs(A (k:n,k))==max(abs(A (k:n,k)))); %寻找主元 t=t(1)+k-1; flag(k)=t; if t~=k
p=A (k,:); A (k,:)=A (t,:); A (t,:)=p; %交换行 det=-det; end if A (k,k)==0
error('矩阵不可逆!'); end
det=det*A (k,k); h=1/A (k,k); A (k,k)=h;
A ([1:k-1 k+1:n],k)=A ([1:k-1 k+1:n],k)*(-h); for i=1:n if i~=k
A (i,[1:k-1 k+1:n])=A (i,[1:k-1 k+1:n])+A (k,[1:k-1…. k+1:n])*A (i,k); %消去 end end
A (k,[1:k-1 k+1:n])=A (k,[1:k-1 k+1:n])*h; end
for k=n-1:-1:1 t=flag(k);
if k~=t %调整 A 列(因为原矩阵换行) p=A (:,t); A (:,t)=A (:,k); A (:,k)=p; end end
8
x18??0.642?x15?例1: 解方程组 ?0.347??0.846x18??0.3x24?751.4x28?230.4x275?90x3.8?4680x3.4?7591x3.2?1?470.4127.7321 1 并验证.
0.8621 解:
>> A=[0.6428 0.3475 -0.8468;0.3475 1.8423 0.4759;-0.8468 0.4759 1.2147]; >> b=[0.4127;1.7321;-0.8621]; >> x=GaJo_inv(A)*b >> A*x-b
(2) Gauss 消去法变形----选主元的三角变形法
通过交换矩阵的行实现矩阵的 PA=LU 分解,采用与列主元相似的方法, 避免一些奇异情况的情况下分解无法进行.具体实现如下: function [L,A,P]=LU_P(A) %选主元的三角变形法P181 %[L,A,P]=LU_P(A) % A 待分解阵 % L 返回单位下三角阵
% A 返回上三角阵,存储在原矩阵中 % P 返回排列阵 n=size(A); if n(1)~=n(2) error('不是方阵!'); end
n=n(1);flag=1:n; P=eye(n,n);L=P; for k=1:n if k~=1
A (k:n,k)=A (k:n,k)-A (k:n,1:k-1)*A (1:k-1,k); end
t=find(abs(A (k:n,k))==max(abs(A (k:n,k)))); t=t(1)+k-1; %选主元 flag(k)=t; if t~=k
p=A (k,:); A (k,:)=A (t,:); A (t,:)=p; end %交换行
9
if abs(A (k,k)) error('Divided by zero!'); end if k~=n A (k+1:n,k)=(1/A (k,k))*A (k+1:n,k); A (k+1:n,k)=A (k+1:n,k); A (k,k+1:n)=A (k,k+1:n)-A (k,1:k-1)*A (1:k-1,k+1:n); end end for k=1:n-1 L (k+1:n,k)=A (k+1:n,k); A (k+1:n,k)=zeros(n-k,1); end %得到 L,U for k=n-1:-1:1 t=flag(k); if k~=t %按主元配制排列阵 p=P(:,t); P(:,t)=P(:,k); P(:,k)=p; end end 例2: 产生 4,4 随机矩阵,做PLU 分解,并作验证. >> A=rand(4,4)*10; >> [L,U,P]=LU_P(A) >> P*A-L*U 数值分析实验二(下) 第三章 解线性方程组的迭代方法 1. 实验目的 用Matlab实现Jacobi 及超松弛跌代法 2. 实验内容 (1) Jacobi 迭代法解线性方程组 10