1.数值模拟MATLAB程序设计
1.1微分方程组模拟、
在微分方程数学模型中，往往需要用到数值模拟。一方面是因为考察模型未来一段时间可能会出现的衰减、震荡或混沌等状况；另一方面，有些微分方程组难以求出解析解，只能求出数值解。
程序代码:
clc;
clear all;
clf;
close all;
% 求解微分方程组
[tt,y] = ode45('lbwfun',0:10:2020,[95.9 0.83 0.003 0.0001 0.0 0.0]);
yx = y(end-10:end,1)
yp = y(end-10:end,2)
yh = y(end-10:end,3)
yt = y(end-10:end,4)
yo = y(end-10:end,5)
ye = y(end-10:end,6)
function ydot = lbwfun(t,y)
%UNTITLED3 此处显示有关此函数的摘要
% 此处显示详细说明
ydot = zeros(6,1);
ydot(1) = 121.793*pi*cos(2*pi*t);
ydot(2) = y(1)-4.03*y(2);
ydot(3) = 0.49*y(2)-17.87*y(3);
ydot(4) = 4.85*y(3)-4.65*y(4);
ydot(5) = 2.55*y(2)+6.12*y(3)+1.95*y(4);
ydot(6) = 1.10*y(2) +6.90*y(3)+2.70*y(4);
end
这里不知道为什么如果把ydot写成矩阵的形式，程序无法执行。就只有分开写了
程序结果：
yx =
95.9240
95.9244
95.9244
95.9247
95.9248
95.9250
95.9251
95.9249
95.9250
95.9252
95.9257
yp =
16.9355
16.9355
16.9356
16.9356
16.9356
16.9356
16.9356
16.9357
16.9357
16.9358
16.9358
yh =
0.4473
0.4473
0.4473
0.4474
0.4473
0.4470
0.4474
0.4473
0.4473
0.4473
0.4473
yt =
0.5807
0.5808
0.5807
0.5807
0.5808
0.5808
0.5807
0.5807
0.5807
0.5807
0.5807
yo =
1.0e+05 *
1.2673
1.2739
1.2805
1.2871
1.2937
1.3003
1.3069
1.3135
1.3201
1.3267
1.3333
ye =
1.0e+04 *
6.2433
6.2758
6.3083
6.3408
6.3734
6.4059
6.4384
6.4709
6.5035
6.5360
6.5685
>>
1.2服从概率分布的随机模拟（排队论）
clear;clc;
% 排队论
% 顾客到达收银台的时间间隔服从平均时间为10s的负指数分布
% 每个顾客的服务时间服从均值为6.5s、标准差为1.2s的正态分布
p = zeros(10,100);
avert = zeros(10,100);
num = 100;
% 分别在顾客人数为10,100,500等情况是，模拟系统工作强度和顾客平均逗留时间
nn = [10 100 500 1000 5000 10000 20000 50000 100000 500000];
for d = 1:length(nn)
for s =1:num
n = nn(d);
dt = exprnd(10,1,n); %到达时间间隔.
st = normrnd(6.5,1.2,1,n); %服务时间间隔
a = zeros(1,n); %每个顾客到达时间
b = zeros(1,n); %开始服务时间
c = zeros(1,n); %结束服务时间
a(1) = 0;
for i = 2:n
a(i) = a(i-1) + dt(i-1);
end
b(1) = 0;
c(1) = st(1);
for i = 2:n
if a(i) > c(i-1)
b(i) = a(i);
else
b(i) = c(i-1);
end
c(i) = b(i) + st(i);
end
% 递推完成，进行统计
p(d,s) = sum(st)/c(end);
avert(d,s) = sum(c-a)/n;
end
end
P = zeros(length(nn),1);
A = P;
for i = 1:length(P)
P(i) = sum(p(i,:))/num;
A(i) = sum(avert(i,:))/num;
end
figure(1);
bar(P);
legend('服务强度');
figure(2);
bar(A);
legend('平均逗留时间');
这里假设顾客不会自动寻优，即不会自动转向空闲的服务台
1.3 蒙特卡洛模拟
1. 蒙特卡洛计算步骤
- 生成随机数
- 利用约束条件检验随机数是否合乎要求
- 求解
2. 求解π
clear;clc;close all;
tic
x = 0:0.01:1;
y = sqrt(1-x.^2);
plot(x,y);
axis square;
grid on;
box on;
format long;
n = 10000000;
count = 0;
for i = 1:n.
point = rand(1,2);
if point(1)^2 + point(2)^2 <=1
count = count+1;
end
end
disp('pi的估计值是:');
count/n*4
toc
3.求解多约束非线性规划
% 利用蒙特卡洛方法求解非线性规划问题
clear;clc;close all;
%产生足够的随机数
N = 10000;
x2 = unifrnd(10,20,N,1);
x3 = unifrnd(-5,16,N,1);
vmax = -inf;
x = zeros(3,1);
for i = 1:N
for j = 1:N
if x2(i)+2*x3(j) >= 10 &&...
3*x2(i)+2*x3(j) <=62;
v = (x2(i)+10)*x2(i)*x3(j);
if v>vmax
vmax = v;
x(1) = x2(i)+10;
x(2) = x2(i);
x(3) = x3(j);
end
end
end
end
x
vmax
2. 动态仿真
1.音频处理
- sound 将向量以声音形式播放(幅度不能超过[-1,1]
- soundsc 将向量整合到[-1,1]区间然后播放
- wavplay 播放.wave格式的音频
- wavrecord 利用Windows音频输入设备进行录音
- wavread 读取.wave格式的音频
2. 动画实现
1.旋转动画
使用rotate命令和drawnow命令
例子:
x = 3*pi*(-1:0.05:1);
y = x;
[X,Y] = meshgrid(x,y);
R = sqrt(X.^2+Y.^2)+eps;
Z = sin(R)./R;
figure(1);
h = surf(X,Y,Z);
axis([-15 15 -15 15 -0.3 1]);
for i= 1:0.0001:2*pi
rotate(h,[0,0,1],i);
drawnow;
pause(0.01);
end
- 彗星轨迹图
- 慧长p*length(y)的二维彗星轨迹：comet(x,y,p)，其中，p的默认值为0.1
- 慧长p*length(y)的三围彗星轨迹:comet3(x,y,z,p)，其中，p的默认值为0.1
3.电影动画
- 调用函数getframe生成每一个帧。
- 调用函数moive按照指定的速度和次数运行该动画
4.实时动画
这是matlab最重要的一种动画。设计实时动画，核心函数指令是erasemode和drawnow。
erasemode指定有4种属性选择：
- normal：重绘整个县市区，它产生的图形准确但是速度慢
- none：不做任何擦除，直接在原图上进行绘画
- xor：擦除旧对象的点并绘制新对象的点
- background：把就对象的颜色变为背景色
指令drawnow使matlab暂停目前的任务而去撒互信屏幕，若不使用指令drawnow，则要等到任务序列执行完毕以后才会刷新屏幕。
这里给出一个小球沿着螺旋线上升的例子：
t = -100:0.1:100;
x = sin(t);
y = cos(t);
z = t;
h = plot3(x,y,z);
% 绘制小球的初始位置
ball = line(x(1),y(1),z(1),'color','red','marker','.','markersize',60);
for i = 1:length(t)
set(ball,'xdata',x(i),'ydata',y(i),'zdata',z(i));
drawnow;
end
网友评论