浙江大学远程教育学院
《可视化计算》课程作业2015年（冬）
姓名： 袁磊 学 号： 715003012008 年级：
15秋计算机
学习中心：
宁波
————————————————————————————— 注意：所有图像的标题必须呈现足够你本人信息
1.（10分）求解下列线性方程组的解： 16
20908356215.87320332231074
445.06.337925.1=-++=++++=++-+-=--+=++++z y v u z y x v u z y x v u z y x v z y x v u
2.（10分）信号y = 5*sin(pi*20*t)+3*cos(2*pi*50*t)幅度为1的白噪声的干扰，请画出此信号，然后进行傅立叶变换，画出变换后的频域信号。
clear
t=0:0.001:0.6;
y=5*sin(pi*20*t)+3*cos(2*pi*50*t);
y=y+randn(1,length(t));
subplot(2,2,1);
plot(y(1:100));
title('袁磊的傅立叶变换频域信号输出图象')
xlabel('变换前的信号');
Y=fft(y,256); subplot(2,2,2);
Y=real(Y); plot(Y(1:256));
xlabel('变换后频域信号')
Pyy=Y.*conj(Y)/256;
f=1000*(0:255)/256;
subplot(2,2,3);
plot(Pyy(1:256));
xlabel('信号功率密度');
3.（10分）在空间有一个带正电的点电荷，请画出此点电荷的空间电位分布与
电场强度的空间分布图。
fprintf('请输入电位分布方程V（x, y)\n');
>>请输入电位分布方程V（x, y)
V= input(':', 's');
:log(x.^2 + y.^2)
NGrid = 20;
xMax = 5;
yMax = 5;
xPlot = linspace(-xMax, xMax, NGrid);
[x, y] = meshgrid(xPlot);
VPlot = eval(V);
[ExPlot, EyPlot] = gradient(-VPlot);
clf;
subplot(1, 2, 1), meshc(VPlot);
set(gcf, 'color', 'w')
xlabel('x');
ylabel('y');
zlabel('电位');
subplot(1, 2, 2), axis([-xMax xMax -yMax yMax]);
cs =contour(x, y, VPlot);
clabel(cs); hold on;
quiver(x, y,ExPlot, EyPlot);
xlabel('x'); ylabel('y'); hold off;
title('电位分布与电位强度线图袁磊制作')
4.（10分）仿照课本第11章的太阳｜地球｜月亮｜卫星，绕转动画实例，呈现地球绕太阳运转的动画。
程序源码：
clear; clc; close all;
x0 = 0; y0 =0; r0=80; Lmin0 = 0; Lmax0 = 0; T0= 2160; w0 = 0*pi/T0; q0 = 0;
x1 = 0; y1 = 0; r1 = 40; Lmin1 = 25; Lmax1 = 30; T1= 1080; w1 = pi/T1; q1 = 0;
hh = figure('numbertitle', 'off','name', '地球绕着太阳转的演示动画');
sun = line(0, 0, 'color', 'r', 'linestyle', '.', 'erasemode', 'xor', 'markersize', r0);
earth = line(x0, y0, 'color', 'k', 'linestyle', '.', 'erasemode', 'xor', 'markersize', r1);
axis off
title('地球绕着太阳转袁磊绘', 'fontname','宋体', 'fontsize', 9, 'FontWeight', 'demi', 'Color', 'black');
text(-50, 50,'太阳');
line(-55, 50, 'color', 'r', 'marker', '.', 'markersize', 80);
text(-50, 40, '地球');
line(-55, 40, 'color', 'k', 'marker', '.', 'markersize', 40);
s1 = [0:.01:2*pi];
line(Lmax1*cos(s1), Lmin1*sin(s1), 'linestyle', ':');
axis([-60, 60, -60, 60]);
t=0;
while 1
if ~ishandle(hh)
return;
end
q0 = t*w0; q1 = t*w1;
t = t+1;
if t>= 4320;
t = 0;
end
x0 = Lmax0 *cos(q1); y1 = Lmin0*sin(q1);
x1 = x0 + Lmax1*cos(q1); y1 = y0 + Lmin1 * sin(q1);
set(sun, 'xdata', x0, 'ydata', y0);
set(earth, 'xdata', x1', 'ydata', y1);
drawnow;
end
5.（10分）设计一个低通滤波器，从混合信号：
x(t)=sin(2*pi*10*t) + cos(2*pi*100*t) + 0.2*randn(size(t)) 中获取10Hz的信号（10分）。
clear;
ws=1000;
t=0:1/ws:0.4;
x=sin(2*pi*10*t)+cos(2*pi*100*t)+0.2*randn(size(t));
wn=ws/2;
[B,A]=butter(10,30/wn);
y=filter(B,A,x);
plot(t,x,'b-')
title('一款低通滤波器--李红辉设计')
hold on
plot(t,y,'r.','MarkerSize',10)
legend('输入','输出')
6.（20分）设计一个程序，应用函数subplot(1,2,1)、subplot(1,2,2)分别显示您本人的二张照片，然后对二张照片分别进行傅立叶变换，并分别画出变换后的频域信号。再把2个频域信号相加，经傅立叶逆变换后，显示时域信号的图像。
A = imread('d:\l1.jpg'); %读取相片图片
B = imread('d:\l2.jpg'); %读取相片图片
subplot(1,2,1); imshow(A)
title('本人第一张照片')
subplot(1,2,2); imshow(B)
title('本人第二张照片')
figure
title('本人二张照片变换后图像')
A1 = fftshift(fft2(A));
B1 = fftshift(fft2(B));
subplot(2,2,1)
Y1 = real(A1);
plot(Y1(1:40))
xlabel('第一张照片变换后频域信号')
subplot(2,2,2)
Y2 = real(B1);
plot(Y2(1:40))
xlabel('第二张照片变换后频域信号')
subplot(2,2,3)
Y3 = Y1 + Y2;
C= ifft2(Y3);
imshow(C);
xlabel('二张照片频域信号相加后图像') （插入相片处）
7.（30分）小论文
根据工作中的实际需要，请设计一个实际工程问题的可视化。可以选择以下之一：（1）工程动画的可视化；（2）大数据处理中的可视化；（3）算法与模型计算的可视化；（4）实际生产流程的可视化；（5）或其它有创新意义的可视化科学计算。要求：
（1）题目有实际意义。
（2）有分析、算法描述
（3）程序源代码设计。
（4）问题结果有可视化显示。
（4）题目的问题有一定的新意。
小论文的字数不能少于2500字，格式由下列各部分组成：
中文题目
摘要：
中文关键词：
英文题目：
英文摘要：
英文关键词：
1.引言
2.算法基础
3.程序代码
4.结果分析（必须有可视化图）
5.结论
参考文献
MATLAB在抛体运动中的探讨
[摘要]计算机在物理教育中的应用已有二十多年的历史，MATLAB 语言是一种集数值计算、
符号运算、可视化建模、仿真和图形处理等多种功能的高级语言。使用MATLAB 模拟物理现象为我们解决问题提供了一种新的方法，利用其方便的数值计算和作图功能，可以方便的模拟一些物理过程。对于处理非线性问题，既能进行数值求解，又能绘制有关曲线，方便实用，基于其功能强大，界面友善，语言自然，交互性强等优点，已成为教学和科研中最基础的软件之一，利用其解决复杂的数值计算问题，可以减少工作量，节约时间，图形绘制问题，真实直观，可以加深理解，提高工作效率。
[关键词]MATLAB 语言 抛体运动 空气阻力 力学 图形绘制
一、 问题的提出
MATLAB 自推向市场以来，得到了广泛的应用和发展，在各高等院校中已经成为线性
代数、自动控制理论、数字信号处理、时间序列分析、动态系统仿真、图像处理等诸多课程的基本教学工具，成为大学生、硕士生、博士生必须掌握的基本技能，尤其在自动控制理论，是最具影响力、最有活力的软件。MATLAB 提供了强大的科学运算、灵活的程序设计流程、高质量的图形可视化与界面设计、便捷的与其它程序和语言接口的功能。对于抛体运动问题，通常需要联立方程组，以及模拟它的路径，运动过程中不同的时间对应不同的位置，利用数学去计算很繁琐，手工绘图误差大，利用MATLAB 可以很好地解决数值计算，模拟抛体运动的路径。
二、 抛体运动的介绍
抛体运动：将物体以一定的初速度向空中抛出，仅在重力作用下物体所作的运动，它的初速度不为零，可分为平抛运动和斜抛运动。物理上提出的“抛体运动”是一种理想化的模型，即把物体看成质点，抛出后只考虑重力作用，忽略空气阻力。抛体运动加速度恒为重力加速度，相等的时间内速度变化量相等，并且速度变化的方向始终是竖直向下的。
一般的处理方法是将其分解为水平方向和竖直方向，平抛运动水平方向是匀速直线运动，竖直方向是自由落体运动，斜抛运动水平方向是匀速直线运动，竖直方向是竖直上抛运动，在任意方向上分解有正交分解和非正交分解两种情加速度及位移等进行相应分析。无论怎样分解，都必须把运动的独立性和独立作用原理结合进行系统分解，即将初速度、受力情、加速度及位移等进行相应分析。
斜抛运动：
水平方向速度
αcos 0
v v x
=
（1）
竖直方向速度
gt v v y
-=αsin 0
（2）
水平方向位移 t x v αcos 0
=
（3）
竖直方向位移 2
021cos gt
t y v -
=
α （4）
平抛运动：
水平方向速度v v x
= （5）
竖直方向速度
gt
v
y
= (6)
水平方向位移t x v 0
=
(7)
竖直方向位移
2
21gt
v y =
(8)
合速度t g v
v v
v y x
t 4220
2
241+
=+=
(9)
合速度方向与水平夹角β：v v
v
gt
tg x
y
==β
(10)
合位移y
x s 2
2
+=
(11)
位移方向与水平夹角α：0
2v gt
tg s
s
x
y
==α
(12)
三、抛体运动的分析
1、斜抛运动的理论分析
忽略空气阻力情况下的抛射体运动是普通物理学中的一个常见问题，在高中物理教材中已有涉及，解决该问题的方法较多，分析的角度也有不同，运用的数学方法也是从初等数学到高等数学而不断深入，对该问题的分析往往是通过运用各种力学原理，推导出该运动的射程，飞行高度，飞行时间以及飞行路径曲线形状等公式，但其数值求解过程比较复杂，因此，对不具备较好高等数学基础的同学来说是较困难的。
设某一抛射体的初速度为0v ，抛射角为θ，将其运动在X,Y 轴上进行正交分解，水平方向速度0cos x v v θ
=
(13) 竖直方向0sin y v v gt
θ=-
(14)
质点的坐标(,)x y 是0()cos()x t t
v θ=
(15)
2
01
()sin 2y t t gt v θ=-
(16)
从上两式消去t ，便得质点的轨迹运动方程2
2
20tan 2cos gx y x v θθ=-
t
(17)
抛射体能达到的最大高度为220
sin 2H g
v
θ=
(18)
其到达最大高度所需时间为0
sin T g
v θ=
(19)
空中飞行时间为 0
sin 22t T g
v
θ
==
(20)
抛射体的最大射程为20
sin 2X g
v θ=
(21)
它跟初速度0v 和抛射角θ有关，在抛射角θ不变的情况下，射程x 与2
0v 成正比，所以射程随初速度的增大而增大。在初速度0v 不变的情况下，随着抛射角θ的增大，射程也增大，当45θ=度时，sin 21θ=，射程达到最大值，以后随着抛射角的增大，射程减小。
利用MATLAB 的绘图功能，可以更直观的体现上述结论。
x=linspace(0,pi/2,100); %产生行向量发射角 g=10; %重力加速度 v1=10; %初速度取10 v2=15;
v3=20; %初速度取20 v4=25; %初速度取25
y1=v1^2*sin(2*x)/g; %初速度为10下的射程 y2=v2^2*sin(2*x)/g; %初速度为15下的射程 y3=v3^2*sin(2*x)/g; %初速度为20下的射程 y4=v4^2*sin(2*x)/g; %初速度为25下的射程 subplot(2,2,1); %选择2*2个区的一号区
plot(x,y1); %输出初速度为10下的射程曲线
title('v0=10'); %加图形标题
text(pi/4,10,'射程为10'); %在最大射程处加图形说明
subplot(2,2,2); %选择2*2个区的二号区
plot(x,y2); %输出初速度为15下的射程曲线
title('v0=15'); %加图形标题
text(pi/4,22.5,'射程为22.5'); %在最大射程处加图形说明
subplot(2,2,3); %选择2*2个区的三号区
plot(x,y3); %输出初速度为20下的射程曲线
title('v0=20'); %加图形标题
text(pi/4,40,'射程为40'); %在最大射程处加图形说明
subplot(2,2,4); %选择2*2个区的四号区
plot(x,y4); %输出初速度为25下的射程曲线
title('v0=25'); %加图形标题
text(pi/4,62.5,'射程为62.5'); %在最大射程处加图形说明
程序运行结果如图1所示。
2、斜抛运动解决实际问题
求解最大飞行路径所对应的抛射角问题（空气阻力忽略不计），如图2所示，
x,处，设在某一微小时段内抛射X,Y坐标轴分别代表抛射体的射程与射高，在()y
体的路径变量
图1 射程与抛射角、初速度的关系
为dt，其对应的水平及竖直方向的变量为dx与dy，
则
22dy dx dL +=
（22）
设射程为
R ，则飞行路径长度
?+=R
dx dx
dy L 0
2
)(
1 （23）
根
据前面的推论，
)
2sin(2
0θg
R v =
(24)
其中v 0为抛射的初始速度，θ为抛射角，
根据运动学原理，有 t x v )cos (0θ= （25）
t gt y v )sin (2
1
02θ+-=
（26） 从（24）、(25）中消除t ，我们可得到该运动的抛物线方程：
θθxtg x g
y v +-
=22
0)
cos (21
（27）
从（24）中可知，为求解L,先得求出
dx
dy
，因此在（4）式两边同时对x 求导，得： θ
θxtg x g
y v +-
=2
0)cos (
（28）
将（27）代入式（24），等式两边同时积分，便得到了飞行路径长度与抛射角之
间的关系：??
?
?????? ??++=θθθθcos sin 1ln cos sin )(22
0g L v
（29）
根据式（28），为求得L 的最大值，将（28）两边同时对θ求导
?
??
??
???? ??+-=θθθθθsin cos 1ln sin 1cos 2)(2
0'g L v
（30）
令0)('=θL ，可得到最大飞行路径所对应的抛射角的大小，但解此方程是比较困
难的。为此，我们采用MATLAB 的函数运算功能来解决这一问题。 程序如下，设其中的抛射初速度s m v 100=，28.9s
m g =。
x=(0:pi/100:pi/2); %产生行向量x
y1=(sin(x)+(cos(x).*cos(x)).*log(1+sin(x))./cos(x))*100/9.8; %飞行路径长度与抛射角之间的函数关系
y2=cos(x).*(1-sin(x).*log((1+sin(x))./cos(x)))*200/9.8;
%飞行路径对抛射角的一阶导数的函数关系
m=(sin(pi/6)+(cos(pi/6)*cos(pi/6))*log(1+sin(pi/6))/cos(pi/6))*100/9.8;
%抛射角取某一特定值时飞行路径值
n=cos(pi/3)*(1-sin(pi/3)*log((1+sin(pi/3))/cos(pi/3)))*200/9.8; %抛射角取某一特定值时飞行路径一阶导的值 plot(x,y1,'b:'); %输出飞行路径长度与抛射角之间的函数表达式 hold on; %设置图形保持状态
plot(x,y2,'k'); % 输出飞行路径对抛射角的一阶导数的函数表达系 hold off; %关闭图形保持
text(pi/6,m,'y1'); %在指定位置添加图例说明 text(pi/3,n,'y2'); %在指定位置添加图列说明 grid; %网格线控制
运行结果如图2所示。
图2给出了飞行路径随抛射角的变化曲线)(θL 及飞行路径曲线的斜度)('θL ，从图中可以得到，当9855.0=θ（弧度）时，即49.56=θ度时，飞行路径最大，
此
时
g L v 2
21.1≈
（31）
我们知道，在不考虑空气阻力的情况下，当抛射角45=θ度时，其射程最远，但此时其飞行路径并不是最远，而是当抛射角49.56=θ度时，其飞行路径最远，且其长度约为g L v 2
021.1≈，实际上，由于空气阻力的存在，抛射体在空中是沿
导弹曲线（弹头飞行时其重心所经过的路线）飞行的，它与抛物线不同，它的升弧与降弧不对称，在重力与空气阻力的共同影响下，弹道形成不均等的圆弧，升弧较长而直伸，降弧较短而弯曲.斜抛射出的炮弹的射程和射高都没有按抛体计算得到的值那么大，路线也不是理想曲线。
图2 抛射角与飞行路径及其一阶导数曲线图
物体在空气中受到的阻力，与物体运动速度大小有密切联系，速度越小，越接近理想情况，当物体速度低于200米每秒时，阻力与物体速度大小的平方成正比，速度介于400至600米每秒之间时，空气阻力与速度大小的三次方成正比，在速度很大的情况下，阻力与速度大小的高次方成正比。 3、抛射角为90度的特殊抛体运动
一弹性小球，初始高度h=10m,向上初速度v0=15米每秒，与地面碰撞的速度衰减系数k=0.8，试计算任意时刻球的位置和速度。
高度与时间的关系：22d y g dt =-，dy
v dt
= （32）
速度与时间关系：
dv
g dt
=-
（33） 对等式两边积分，有dv gdt =-??
，0v v gt =- （34）
2001
2
y y t gt v =+- （35）
由此可得数学方程：
第一次落地前：01v v gt =- （36）
2
012
gt y h v t =+-
（37）
1 3.62T s = （38）
第二次落地前：
02011()v k v gT =-- （39） 02v v gt =-
（40）
2
022
gt y v t =-
(41）
02
22v T g
=
(42)
第三次落地前：
03022()v k v gT =-- （43） 03v v gt =- （44）
2032
gt y v t =-
（45）
03
32v T g
=
（46）
. . . . . .
第n 次落地前：
00(1)(1)()n n n v k v gT --=-- （47）
0n v v gt =- （48）
202
m gt y v =-
（49)
02n
n v T g
=
(50) 如用手工进行计算，计算量极大，利用MATLAB 编程如下：
v0=15; %初速度 h=10; %初始高度 g=-9.8; %重力加速度 k=0.8; %衰减系数 T=0; %落地时间
for t=0:0.05:20 % 产生时间的行向量 v=v0+g*(t-T); %求速度 y=h+v0*(t-T)+g*(t-T)^2/2; %求高度 if y<=0 %循环判断条件 v0=-k*v; %衰减的速度 T=t; %求球每次落地所用时间 h=0; %将高度变零 end %选择结构结束
subplot(1,2,1); %选择1*2中的一号区 pause(0.1); %延缓
plot(1,y,'or','MarkerSize',10,'Markerface',[1,0,0]);
%输出求球的运动图像
title('运动变化图'); %图形名称
axis([0,2,0,25]); %坐标控制
subplot(2,2,2); %选择2*2中的二号区
axis([0,20,-25,30]); %坐标控制
grid on; %不画网格线
plot(t,v,'*r','MarkerSize',2); %画球的速度曲线
xlabel('时间t'); %坐标轴说明
ylabel('速度v'); %坐标轴说明
title('速度变化趋势图'); %图形名称
hold on; %设置图形保持状态
subplot(2,2,4); %选择2*2中的四号区
axis([0,20,0,25]); %坐标控制
grid on; %不加网格线
plot(t,y,'*b','MarkerSize',2); %画球的位置曲线
xlabel('时间t'); %坐标轴说明
ylabel('高度y'); %坐标轴说明
title('位置变化图'); %图形名称
grid on %不加网格线
hold on %设置图形保持状态
end %循环结束
程序运行结果如图3所示。
4、对平抛运动的分析
将物体用一定的初速度沿水平方向抛出，不考虑空气的阻力，物体只在重力作用下所做的运动，叫做平抛运动。
竖直的重力与速度方向有夹角，做曲线运动；水平方向不受外力作用，是匀速运动，速度为Vo；竖直方向受重力作用，没有初速度，加速度为重力加速度g,是自由体运动。即做平抛运动的物体，在水平方向上由于不受力，将作匀速直线运动；在竖直方向上的物体的初速度为0.且只受到重力作用，物体做自由落体运动，加速度为g。
平抛运动的规律：
1、抛出t秒末的速度：
一抛出点为坐标原点，水平方向为x轴（正方向和初速度V0的方向相同），
竖直方向为y轴，正方向向下，则：
水平速度：Vx=Vo （51）竖直分速度：Vy=gt （52）
合速度：Vt=（53）
tanθ=Vy
Vx
=
gt
Vo
（54）
运用MATLAB编程得到速度随时间的变化关系，程序如下：
图3 小球落地速度及位置曲线
t=0:0.01:10; %产生时间的行向量 Vt=-sqrt(10^2+9.8*t.^2); %求速度 plot(t,Vt); %输出速度曲线 title('物体速度随时间的变化'); % 图形名称 grid %加网格线
运行结果如图4所示。
2、平抛运动的物体在任意时刻t 的位置坐标：
水平位移：x=Vot （55） 竖直位移：y=
12
g 2
t （56） 合位移：2
2
x y + （57）
tan θ=
y x =2gt
Vo （58） 运用该公式，我们可以求得物体在任意时刻的坐标并找到物体所在位置后，再用平滑曲线把这些点连起来，就得到平抛运动的轨迹。
运用MATLAB 编程的到物体运动的曲线的程序如下： t=0:0.01:10; %产生时间行向量
s=-sqrt((3*t).^2+(0.5*9.8*t.^2).^2); %求位移 plot(t,s,'r:'); %输出位移曲线 title('物体平抛运动轨迹'); %图形名称
grid %加网格线运行结果如图5所示。
图4 平抛运动速度随时间变化关系
图5 物体平抛轨迹曲线
四、结论
从以上对抛体运动的分析可得出这些结论：
1、抛射体的射程与初速度和抛射角有关，在抛射角不变的情况下，射程随初速度的增大而增大，在抛射角不变的情况下，射程随抛射角的增大而增大，当抛射角达到四十五度时射程达到最大值，之后射程随着抛射角的增大而减小。
2、抛体运动的分析牵扯到很多复杂的推理与计算，MATLAB可以有效的解决这一问题，
MATLAB的数值计算和作图功能来模拟物理现象，减少了计算工作量，并且可以准确的得到计算结果，且绘制的图形使抽象的问题形象化。
MATLAB具有程序控制结构、函数调用、数据结构、输入输出、面向对象等程序语言特征，因此与传统编程语言一样，可以进行程序设计，而且简单易学，编程效率高。MATLAB以矩阵最为数据操作的基本单位，使得矩阵运算变得简捷、方便、高效，其还提供了十分丰富的数值计算函数、符号计算函数。用其绘图十分方便，还可对图形进行修饰和控制。
[参考文献]
[1] 刘卫国.MATLAB程序设计与应用 (第二版)[M].北京：高等教育出版社，2006.
[2] 同济大学数学系编.高等数学（第六版）上册[M].北京：高等教育出版社，2007.
[3]马文蔚.物理学（上册）（第四版）[M].北京：高等教育出版社，1999.