根者稿逗据Doolittle分解格式， 可以用matlab分解A矩阵为一个下三角矩阵L与首卖上三角矩阵U的乘积。
其实现过程为：
第一步：初始化
1、初始化上三角阵的第一行
2、初始化下三敬薯角阵的第一列
第二步：前向分解计算
根据上述方法编程，为了说明问题，特举例如下。
3、例题
将 A=[1 11 2]分解成L和U矩阵
将 A =[2,2,34,7,7-2,4,5]分解成L和U矩阵
4、运行程序，可得到如下结果。
1-32是：图形应用篇33-66是：界面设计篇
67-84是：图形处理篇
85-100是：数值分析篇
function shili01
h0=figure('toolbar','none',...
'position',[198 56 350 300],...
'name','实例01')
h1=axes('parent',h0,...
'visible','off')
x=-pi:0.05:pi
y=sin(x)
plot(x,y)
xlabel('自变量X')
ylabel('函数值Y')
title('SIN( )函数曲线')
grid on
实例2：三角函纤梁数曲线（2）
function shili02
h0=figure('toolbar','none',...
'position',[200 150 450 350],...
'name','实例02')
x=-pi:0.05:pi
y=sin(x)+cos(x)
plot(x,y,'-*r','linewidth'陪宽,1)
grid on
xlabel('自变量X')
ylabel('函数值Y')
title('三角函数')
实例3：图形的叠加
function shili03
h0=figure('toolbar','none',...
'position',[200 150 450 350],...
'name','实例03')
x=-pi:0.05:pi
y1=sin(x)
y2=cos(x)
plot(x,y1,...
'-*r',...
x,y2,...
'--og')
grid on
xlabel('自变量X')
ylabel('函数值Y')
title('三角函数毁乱运')
实例4：双y轴图形的绘制
function shili04
h0=figure('toolbar','none',...
'position',[200 150 450 250],...
'name','实例04')
x=0:900a=1000b=0.005
y1=2*x
y2=cos(b*x)
[haxes,hline1,hline2]=plotyy(x,y1,x,y2,'semilogy','plot')
axes(haxes(1))
ylabel('semilog plot')
axes(haxes(2))
ylabel('linear plot')
实例5：单个轴窗口显示多个图形
function shili05
h0=figure('toolbar','none',...
'position',[200 150 450 250],...
'name','实例05')
t=0:pi/10:2*pi
[x,y]=meshgrid(t)
subplot(2,2,1)
plot(sin(t),cos(t))
axis equal
subplot(2,2,2)
z=sin(x)-cos(y)
plot(t,z)
axis([0 2*pi -2 2])
subplot(2,2,3)
h=sin(x)+cos(y)
plot(t,h)
axis([0 2*pi -2 2])
subplot(2,2,4)
g=(sin(x).^2)-(cos(y).^2)
plot(t,g)
axis([0 2*pi -1 1])
实例6：图形标注
function shili06
h0=figure('toolbar','none',...
'position',[200 150 450 400],...
'name','实例06')
t=0:pi/10:2*pi
h=plot(t,sin(t))
xlabel('t=0到2\pi','fontsize',16)
ylabel('sin(t)','fontsize',16)
title('\it{从 0to2\pi 的正弦曲线}','fontsize',16)
x=get(h,'xdata')
y=get(h,'ydata')
imin=find(min(y)==y)
imax=find(max(y)==y)
text(x(imin),y(imin),...
['\leftarrow最小值=',num2str(y(imin))],...
'fontsize',16)
text(x(imax),y(imax),...
['\leftarrow最大值=',num2str(y(imax))],...
'fontsize',16)
实例7：条形图形
function shili07
h0=figure('toolbar','none',...
'position',[200 150 450 350],...
'name','实例07')
tiao1=[562 548 224 545 41 445 745 512]
tiao2=[47 48 57 58 54 52 65 48]
t=0:7
bar(t,tiao1)
xlabel('X轴')
ylabel('TIAO1值')
h1=gca
h2=axes('position',get(h1,'position'))
plot(t,tiao2,'linewidth',3)
set(h2,'yaxislocation','right','color','none','xticklabel',[])
实例8：区域图形
function shili08
h0=figure('toolbar','none',...
'position',[200 150 450 250],...
'name','实例08')
x=91:95
profits1=[88 75 84 93 77]
profits2=[51 64 54 56 68]
profits3=[42 54 34 25 24]
profits4=[26 38 18 15 4]
area(x,profits1,'facecolor',[0.5 0.9 0.6],...
'edgecolor','b',...
'linewidth',3)
hold on
area(x,profits2,'facecolor',[0.9 0.85 0.7],...
'edgecolor','y',...
'linewidth',3)
hold on
area(x,profits3,'facecolor',[0.3 0.6 0.7],...
'edgecolor','r',...
'linewidth',3)
hold on
area(x,profits4,'facecolor',[0.6 0.5 0.9],...
'edgecolor','m',...
'linewidth',3)
hold off
set(gca,'xtick',[91:95])
set(gca,'layer','top')
gtext('\leftarrow第一季度销量')
gtext('\leftarrow第二季度销量')
gtext('\leftarrow第三季度销量')
gtext('\leftarrow第四季度销量')
xlabel('年','fontsize',16)
ylabel('销售量','fontsize',16)
实例9：饼图的绘制
function shili09
h0=figure('toolbar','none',...
'position',[200 150 450 250],...
'name','实例09')
t=[54 21 35
68 54 35
45 25 12
48 68 45
68 54 69]
x=sum(t)
h=pie(x)
textobjs=findobj(h,'type','text')
str1=get(textobjs,{'string'})
val1=get(textobjs,{'extent'})
oldext=cat(1,val1{:})
names={'商品一：''商品二：''商品三：'}
str2=strcat(names,str1)
set(textobjs,{'string'},str2)
val2=get(textobjs,{'extent'})
newext=cat(1,val2{:})
offset=sign(oldext(:,1)).*(newext(:,3)-oldext(:,3))/2
pos=get(textobjs,{'position'})
textpos=cat(1,pos{:})
textpos(:,1)=textpos(:,1)+offset
set(textobjs,{'position'},num2cell(textpos,[3,2]))
实例10：阶梯图
function shili10
h0=figure('toolbar','none',...
'position',[200 150 450 400],...
'name','实例10')
a=0.01
b=0.5
t=0:10
f=exp(-a*t).*sin(b*t)
stairs(t,f)
hold on
plot(t,f,':*')
hold off
glabel='函数e^{-(\alpha*t)}sin\beta*t的阶梯图'
gtext(glabel,'fontsize',16)
xlabel('t=0:10','fontsize',16)
axis([0 10 -1.2 1.2])
实例11：枝干图
function shili11
h0=figure('toolbar','none',...
'position',[200 150 450 350],...
'name','实例11')
x=0:pi/20:2*pi
y1=sin(x)
y2=cos(x)
h1=stem(x,y1+y2)
hold on
h2=plot(x,y1,'^r',x,y2,'*g')
hold off
h3=[h1(1)h2]
legend(h3,'y1+y2','y1=sin(x)','y2=cos(x)')
xlabel('自变量X')
ylabel('函数值Y')
title('正弦函数与余弦函数的线性组合')
实例12：罗盘图
function shili12
h0=figure('toolbar','none',...
'position',[200 150 450 250],...
'name','实例12')
winddirection=[54 24 65 84
256 12 235 62
125 324 34 254]
windpower=[2 5 5 3
6 8 12 7
6 14 10 8]
rdirection=winddirection*pi/180
[x,y]=pol2cart(rdirection,windpower)
compass(x,y)
desc={'风向和风力',
'北京气象台',
'10月1日0：00到',
'10月1日12：00'}
gtext(desc)
实例13：轮廓图
function shili13
h0=figure('toolbar','none',...
'position',[200 150 450 250],...
'name','实例13')
[th,r]=meshgrid((0:10:360)*pi/180,0:0.05:1)
[x,y]=pol2cart(th,r)
z=x+i*y
f=(z.^4-1).^(0.25)
contour(x,y,abs(f),20)
axis equal
xlabel('实部','fontsize',16)
ylabel('虚部','fontsize',16)
h=polar([0 2*pi],[0 1])
delete(h)
hold on
contour(x,y,abs(f),20)
实例14：交互式图形
function shili14
h0=figure('toolbar','none',...
'position',[200 150 450 250],...
'name','实例14')
axis([0 10 0 10])
hold on
x=[]
y=[]
n=0
disp('单击鼠标左键点取需要的点')
disp('单击鼠标右键点取最后一个点')
but=1
while but==1
[xi,yi,but]=ginput(1)
plot(xi,yi,'bo')
n=n+1
disp('单击鼠标左键点取下一个点')
x(n,1)=xi
y(n,1)=yi
end
t=1:n
ts=1:0.1:n
xs=spline(t,x,ts)
ys=spline(t,y,ts)
plot(xs,ys,'r-')
hold off
实例14：交互式图形
function shili14
h0=figure('toolbar','none',...
'position',[200 150 450 250],...
'name','实例14')
axis([0 10 0 10])
hold on
x=[]
y=[]
n=0
disp('单击鼠标左键点取需要的点')
disp('单击鼠标右键点取最后一个点')
but=1
while but==1
[xi,yi,but]=ginput(1)
plot(xi,yi,'bo')
n=n+1
disp('单击鼠标左键点取下一个点')
x(n,1)=xi
y(n,1)=yi
end
t=1:n
ts=1:0.1:n
xs=spline(t,x,ts)
ys=spline(t,y,ts)
plot(xs,ys,'r-')
hold off
实例15：变换的傅立叶函数曲线
function shili15
h0=figure('toolbar','none',...
'position',[200 150 450 250],...
'name','实例15')
axis equal
m=moviein(20,gcf)
set(gca,'nextplot','replacechildren')
h=uicontrol('style','slider','position',...
[100 10 500 20],'min',1,'max',20)
for j=1:20
plot(fft(eye(j+16)))
set(h,'value',j)
m(:,j)=getframe(gcf)
end
clf
axes('position',[0 0 1 1])
movie(m,30)
实例16：劳伦兹非线形方程的无序活动
function shili15
h0=figure('toolbar','none',...
'position',[200 150 450 250],...
'name','实例15')
axis equal
m=moviein(20,gcf)
set(gca,'nextplot','replacechildren')
h=uicontrol('style','slider','position',...
[100 10 500 20],'min',1,'max',20)
for j=1:20
plot(fft(eye(j+16)))
set(h,'value',j)
m(:,j)=getframe(gcf)
end
clf
axes('position',[0 0 1 1])
movie(m,30)
实例17：填充图
function shili17
h0=figure('toolbar','none',...
'position',[200 150 450 250],...
'name','实例17')
t=(1:2:15)*pi/8
x=sin(t)
y=cos(t)
fill(x,y,'r')
axis square off
text(0,0,'STOP',...
'color',[1 1 1],...
'fontsize',50,...
'horizontalalignment','center')
例18：条形图和阶梯形图
function shili18
h0=figure('toolbar','none',...
'position',[200 150 450 250],...
'name','实例18')
subplot(2,2,1)
x=-3:0.2:3
y=exp(-x.*x)
bar(x,y)
title('2-D Bar Chart')
subplot(2,2,2)
x=-3:0.2:3
y=exp(-x.*x)
bar3(x,y,'r')
title('3-D Bar Chart')
subplot(2,2,3)
x=-3:0.2:3
y=exp(-x.*x)
stairs(x,y)
title('Stair Chart')
subplot(2,2,4)
x=-3:0.2:3
y=exp(-x.*x)
barh(x,y)
title('Horizontal Bar Chart')
实例19：三维曲线图
function shili19
h0=figure('toolbar','none',...
'position',[200 150 450 400],...
'name','实例19')
subplot(2,1,1)
x=linspace(0,2*pi)
y1=sin(x)
y2=cos(x)
y3=sin(x)+cos(x)
z1=zeros(size(x))
z2=0.5*z1
z3=z1
plot3(x,y1,z1,x,y2,z2,x,y3,z3)
grid on
xlabel('X轴')
ylabel('Y轴')
zlabel('Z轴')
title('Figure1:3-D Plot')
subplot(2,1,2)
x=linspace(0,2*pi)
y1=sin(x)
y2=cos(x)
y3=sin(x)+cos(x)
z1=zeros(size(x))
z2=0.5*z1
z3=z1
plot3(x,z1,y1,x,z2,y2,x,z3,y3)
grid on
xlabel('X轴')
ylabel('Y轴')
zlabel('Z轴')
title('Figure2:3-D Plot')
实例20：图形的隐藏属性
function shili20
h0=figure('toolbar','none',...
'position',[200 150 450 300],...
'name','实例20')
subplot(1,2,1)
[x,y,z]=sphere(10)
mesh(x,y,z)
axis off
title('Figure1:Opaque')
hidden on
subplot(1,2,2)
[x,y,z]=sphere(10)
mesh(x,y,z)
axis off
title('Figure2:Transparent')
hidden off
实例21PEAKS函数曲线
function shili21
h0=figure('toolbar','none',...
'position',[200 100 450 450],...
'name','实例21')
[x,y,z]=peaks(30)
subplot(2,1,1)
x=x(1,:)
y=y(:,1)
i=find(y>0.8&y<1.2)
j=find(x>-0.6&x<0.5)
z(i,j)=nan*z(i,j)
surfc(x,y,z)
xlabel('X轴')
ylabel('Y轴')
zlabel('Z轴')
title('Figure1:surfc函数形成的曲面')
subplot(2,1,2)
x=x(1,:)
y=y(:,1)
i=find(y>0.8&y<1.2)
j=find(x>-0.6&x<0.5)
z(i,j)=nan*z(i,j)
surfl(x,y,z)
xlabel('X轴')
ylabel('Y轴')
zlabel('Z轴')
title('Figure2:surfl函数形成的曲面')
实例22：片状图
function shili22
h0=figure('toolbar','none',...
'position',[200 150 550 350],...
'name','实例22')
subplot(1,2,1)
x=rand(1,20)
y=rand(1,20)
z=peaks(x,y*pi)
t=delaunay(x,y)
trimesh(t,x,y,z)
hidden off
title('Figure1:Triangular Surface Plot')
subplot(1,2,2)
x=rand(1,20)
y=rand(1,20)
z=peaks(x,y*pi)
t=delaunay(x,y)
trisurf(t,x,y,z)
title('Figure1:Triangular Surface Plot')
实例23：视角的调整
function shili23
h0=figure('toolbar','none',...
'position',[200 150 450 350],...
'name','实例23')
x=-5:0.5:5
[x,y]=meshgrid(x)
r=sqrt(x.^2+y.^2)+eps
z=sin(r)./r
subplot(2,2,1)
surf(x,y,z)
xlabel('X-axis')
ylabel('Y-axis')
zlabel('Z-axis')
title('Figure1')
view(-37.5,30)
subplot(2,2,2)
surf(x,y,z)
xlabel('X-axis')
ylabel('Y-axis')
zlabel('Z-axis')
title('Figure2')
view(-37.5+90,30)
subplot(2,2,3)
surf(x,y,z)
xlabel('X-axis')
ylabel('Y-axis')
zlabel('Z-axis')
title('Figure3')
view(-37.5,60)
subplot(2,2,4)
surf(x,y,z)
xlabel('X-axis')
ylabel('Y-axis')
zlabel('Z-axis')
title('Figure4')
view(180,0)
实例24：向量场的绘制
function shili24
h0=figure('toolbar','none',...
'position',[200 150 450 350],...
'name','实例24')
subplot(2,2,1)
z=peaks
ribbon(z)
title('Figure1')
subplot(2,2,2)
[x,y,z]=peaks(15)
[dx,dy]=gradient(z,0.5,0.5)
contour(x,y,z,10)
hold on
quiver(x,y,dx,dy)
hold off
title('Figure2')
subplot(2,2,3)
[x,y,z]=peaks(15)
[nx,ny,nz]=surfnorm(x,y,z)
surf(x,y,z)
hold on
quiver3(x,y,z,nx,ny,nz)
hold off
title('Figure3')
subplot(2,2,4)
x=rand(3,5)
y=rand(3,5)
z=rand(3,5)
c=rand(3,5)
fill3(x,y,z,c)
grid on
title('Figure4')
实例25：灯光定位
function shili25
h0=figure('toolbar','none',...
'position',[200 150 450 250],...
'name','实例25')
vert=[1 1 11 2 1
2 2 12 1 1
1 1 21 2 2
2 2 22 1 2]
fac=[1 2 3 42 6 7 3
4 3 7 81 5 8 4
1 2 6 55 6 7 8]
grid off
sphere(36)
h=findobj('type','surface')
set(h,'facelighting','phong',...
'facecolor',...
'interp',...
'edgecolor',[0.4 0.4 0.4],...
'backfacelighting',...
'lit')
hold on
patch('faces',fac,'vertices',vert,...
'facecolor','y')
light('position',[1 3 2])
light('position',[-3 -1 3])
material shiny
axis vis3d off
hold off
实例26：柱状图
function shili26
h0=figure('toolbar','none',...
'position',[200 50 450 450],...
'name','实例26')
subplot(2,1,1)
x=[5 2 1
8 7 3
9 8 6
5 5 5
4 3 2]
bar(x)
xlabel('X轴')
ylabel('Y轴')
title('第一子图')
subplot(2,1,2)
y=[5 2 1
8 7 3
9 8 6
5 5 5
4 3 2]
barh(y)
xlabel('X轴')
ylabel('Y轴')
title('第二子图')
1.新建一个m文件，存储被积函数：function f=fun(x)
f=1./x
2.另建一个m文件，放入Romberg积分程序：
function [I,T]=Romberg(fun,a,b,eps,varargin)
% Romberg公式求解数值积分
% 输入参数：
% ---fun：被积函数
% ---a,b：积分区间的端点
% ---eps：精度要派告求，默认值为1e-6
% ---p1,p2,...：fun的附加参数
% 输出参数：
% ---I：求得的积分值
% ---T：Romberg积分过程产生的下三角阵
if nargin<4|isempty(eps)
eps=1e-6
end
N=1h=b-a
T(1,1)=h/2*sum(feval(fun,[a,b],varargin{:}))
tol=1
while tol>eps
h=h/2N=2*Nk=log2(N)
x=a+(2*(1:N/2)-1)*h
fx = feval(fun,x,varargin{:})% 计算函数值
T(k+1,1)=1/2*T(k,1)+h*sum(fx)
for j=1:k
T(k+1,j+1)=(4^j*T(k+1,j)-T(k,j))/(4^j-1)
end
tol=abs(T(k+1,k+1)-T(k,k))
end
I=T(k+1,k+1)
3.在命令窗口调用计算：
>>y=vpa(Romberg(@fun,1,3,1e-5),6)
y =
1.09861
>>烂羡则饥棚
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）