history080501.m

(有源代码)数值分析作业,本文主要包括两个部分,第一部分是常微分方程(ODE)的三个实验题,第二部分是有关的拓展讨论,包括高阶常微分的求解和边值问题的求解(BVP).文中的算法和算例都是基于Matla

function history080501
% history of 080501

%-- Unknown date --%
afyy=fyy.*conj(fyy)/4096;
plot(f,afyy(1:2049))
fyy=fft(yyf,4096);
afyy=fyy.*conj(fyy)/4096;
plot(f,afyy(1:2049))
plot(xx,yyf)
plot(xx,yy)
plot(xx,yy,xx,yyf)
plot(xx,yy,xx,yyf,xx,sin(2*pi*5*xx))
Fs = 8192;  % Sampling Frequency
N   = 10;    % Order
Fc1 = 3;    % First Cutoff Frequency
Fc2 = 1000;  % Second Cutoff Frequency
h  = fdesign.bandpass('N,Fc1,Fc2', N, Fc1, Fc2, Fs);
Hd = butter(h);
yyf=filter(Hd,yy);
hold off
plot(xx,yy,xx,yyf,xx,sin(2*pi*5*xx))
fyy=fft(yyf,4096);
afyy=fyy.*conj(fyy)/4096;
plot(f,afyy(1:2049))
absfyy=abs(fyy);
plot(f,afyy(1:2049))
doc buttord
doc freqz
freqz(Hd)
Fs = 8192;  % Sampling Frequency
N   = 10;    % Order
Fc1 = 3;    % First Cutoff Frequency
Fc2 = 4000;  % Second Cutoff Frequency
h  = fdesign.bandpass('N,Fc1,Fc2', N, Fc1, Fc2, Fs);
Hd = butter(h);
freqz(Hd)
[H,W]=freqz(Hd)
clc
[H,W]=freqz(Hd);
plot (W, 'DisplayName', 'W', 'YDataSource', 'W'); figure(gcf)
plot (H, 'DisplayName', 'H', 'YDataSource', 'H'); figure(gcf)
plot (W, 'DisplayName', 'W', 'YDataSource', 'W'); figure(gcf)
plot(W*8192/(2*pi),abs(H));grid;
plot (W, 'DisplayName', 'W', 'YDataSource', 'W'); figure(gcf)
grid
load dtdat
doc psd
doc signal/fdatool
wvtool(hanning(64))
help rect
wvtool(rect(64))
wvtool(rectwin(64))
wvtool(rectwin(16))
wvtool(rectwin(128))
doc specgram
disp(Hd)
plot (run16t_data, 'DisplayName','run16t_data', 'YDataSource', 'run16t_data'); figure(gcf)
plot (ty, 'DisplayName', 'ty', 'YDataSource', 'ty'); figure(gcf)
plot (x, 'DisplayName', 'x', 'YDataSource', 'x'); figure(gcf)
plot (y, 'DisplayName', 'y', 'YDataSource', 'y'); figure(gcf)
clc
x=run16t_data(:,1);
y=run16t_data(:,2);
Fs=8192;
sptool
load mtlb
cc
load mtlb
sptool
load wjl0313
plot (y, 'DisplayName', 'y', 'YDataSource', 'y'); figure(gcf)
Fs=8192
t = 0:0.001:0.6;
x = sin(2*pi*50*t)+sin(2*pi*120*t);
y = x + 2*randn(size(t));
Fs=1000;
fy=fft(y,512);
f=1000*(0:256)/512;
plot(f,fy(1:257))
plot(f,abs(fy(1:257)))
plot(f,abs(fy(1:257))/256)
pyy=fy.*conj(fy)/256;
plot(f,pyy(1:257))
pyy=fy.*conj(fy)/512;
spyy=sqrt(pyy);
plot(f,spyy(1:257))
doc fft
randn('state',0);
fs = 1000;           % Sampling frequency
t = (0:fs)/fs;       % One second worth of samples
A = [1 2];           % Sinusoid amplitudes (row vector)
f = [150;140];       % Sinusoid frequencies (column vector)
xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));
Hs = spectrum.periodogram('Hamming');
psd(Hs,xn,'Fs',fs,'NFFT',1024,'SpectrumType','twosided')
[Pxx,F] = psd(Hs,xn,fs,'twosided');
Pow = (fs/length(Pxx)) * sum(Pxx)
[Pxx,F] = psd(Hs,xn,fs,'twosided');
[Pxxo,F] = psd(Hs,xn,fs,'onesided');
randn('state',1)
fs = 1000;               % Sampling frequency
t = (0:0.3*fs)./fs;      % 301 samples
A = [2 8];               % Sinusoid amplitudes (row vector)
f = [150;140];           % Sinusoid frequencies (column vector)
xn = A*sin(2*pi*f*t) + 5*randn(size(t));
Hs = spectrum.periodogram('rectangular')
psd(Hs,xn,'Fs',fs,'NFFT',1024);
Hs = spectrum.welch('rectangular',150,50);
psd(Hs,xn,'Fs',fs,'NFFT',512);
randn('state',0)
fs = 1000;               % Sampling frequency
t = (0:fs)/fs;           % One second worth of samples
A = [1 2];               % Sinusoid amplitudes
f = [150;140];           % Sinusoid frequencies
xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));
Hs1 = spectrum.mtm(4,'adapt');
psd(Hs1,xn,'Fs',fs,'NFFT',1024);
figure
plot(t,xn)
[Pxx1,F] = psd(Hs1,xn,fs);
h = ones(1,10)/10;
h
yn = filter(h,1,xn)
clc
[HEST,f] = tfestimate(xn,yn,256,128,256,fs);
H = freqz(h,1,f,fs);
plot (H, 'DisplayName', 'H', 'YDataSource', 'H'); figure(gcf)
subplot(2,1,1); plot(f,abs(H));
title('Actual Transfer Function Magnitude');
subplot(2,1,2); plot(f,abs(HEST));
title('Transfer Function Magnitude Estimate');
xlabel('Frequency (Hz)');
mscohere(xn,yn,256,128,256,fs)
%-- 08-3-15 上午9:35 --%
whats
clc
what
load fun16t_data
clear
load run16t
plot(run16t_data(:,2))
y = run16t_data(:,2);
x = run16t_data(:,1);
plot(x,y)
sptool
yf=run16t_f.data;
plot(x,y,x,yf)
plot(x,y,x,yf,'r')
figure
plot(x,y,x,yf,'r')
figure
plot(x,y-yf);
fye=fft(y-yf);
plot(abs(fye))
figure
plot(abs(fft(y)))
hold on
plot(abs(fye),'r')
plot(abs(fft(yf)),'y')
figure(1)
plot(x,y)
y=9.8 * y;
plot(x,y)
figure
plot(x,y)
hold on
plot(x,yf,'r')
what
yv=trapz(x,y);
figure
plot(x,yv)
plot(x(1:end-1),yv)
ys=trapz(x,yv);
ys=trapz(x(1:end-1),yv);
figure
plot(x(1:end-2),ys);
hold on
plot(x(1:end-1),yv)
plot(x,y)
figure
plot(x,y)
subplot(1,3,1)
plot(x,y)
subplot(1,3,2)
plot(x(1:end-1),yv)
subplot(1,3,3)
plot(x(1:end-2),ys);
figure
subplot(1,3,1)
plot(x,y)
subplot(1,3,2)
plot(x(1:end-1),yv)
subplot(1,3,3)
plot(x(1:end-2),ys);
figure
subplot(3,1,1)
plot(x,y)
subplot(3,1,2)
plot(x(1:end-1),yv)
subplot(3,1,3)
plot(x(1:end-2),ys);
help detrends
help detrend
who
whos
yvdetrend=detrend(yv);
figure
plot(x(1:end-1),yv)
plot(x(1:end-1),yvdetrend)
hold on
plot(x(1:end-1),yv)
hold off
plot(x(1:end-1),yvdetrend,x(1:end-1),yv,'r')
plot(x(1:end-1),yvdetrend,'r',x(1:end-1),yv)
title('yv detrend and yv ')
figure
plot(abs(fft(yv)))
figure(1)
plot(abs(fft(yvdetrend)))
grid
figure(1)
plot(x(1:end-2),ys);
title('ys')
doc polyfit
xx=x(1:end-2);
yy=ys;
plot(xx,yy)
p=polyfit(xx,yy,2)
ff=polyval(p,xx);
ysdetrend=yy-ff;
plot(xx,yy,xx,ff,'r')
figure
plot(xx,yy,xx,ysdetrend,'r')
plot(xx,ysdetrend)
hold on
10e-7
1e-7
plot(xx,1e-7)
plot(xx,1e-6)
plot(xx,-1e-6)
clc
plot(xx,ysdetrend)
plot(xx,1e-6)
plot(xx,-1e-6)
hold on
plot(xx,1e-6)
plot(xx,ysdetrend)
plot(xx,ys)
save wjlyyyvys all
save wjlyyyvys
what
clear
cc
xx=0:0.1:10*pi;
yy=sin(xx);
yyv=trapz(xx,yy);
plot(xx,yy,xx(1:end-1),yyv,'r')
yys=trapz(xx(1:end-1),yyv);
figure
plot(xx,yy,xx(1:end-1),yy,'r',xx(1:end-2),yys,'y'))
plot(xx,yy,xx(1:end-1),yy,'r',xx(1:end-2),yys,'y')
figure
plot(xx,yy,xx(1:end-1),yy,'r',xx(1:end-2),yys,'y')
plot(xx(1:end-2),yys)
plot(xx(1:end-2),yys,'y',xx,yy,xx(1:end-1),yy,'r'))
plot(xx(1:end-2),yys,'y',xx,yy,xx(1:end-1),yy,'r')
plot(xx(1:end-2),yys,'y',xx,yy,'b',xx(1:end-1),yy,'r')
clf
plot(xx(1:end-2),yys)
hold on
plot(xx(1:end-1),yyv)
plot(xx,yy)
xx1=xx(1:end-1);
xx2=xx(1:end-2);
yyvd=detrend(xx1,yyv);
yyvd=detrend(yyv);
plot(xx1,yyvd,'g');
yys1=trapz(xx1,yyvd);
yysd=detrend(yys1);
plot(xx2,yysd,'b')
plot(xx2,yysd,'k')
p=polyfit(xx2,yys,2)
p=polyfit(xx2',yys,2)
pp=polyval(p,xx2');
yysdd=yys-pp;
plot(xx2,yysdd,'-.')
plot(xx2,pp,'-.')
cc
13.85/31.2*1000
443.9/1.4
31.2*1.4
ans/281
ans*2000
cc
xx=0:0.1:10*pi;
yy=sin(2*xx) +rand(xx);
yy=sin(2*xx) +rand(size(xx));
plot(xx,yy)
xx1=xx(1:end-1);
xx2=xx(1:end-2);
yv=trapz(xx,yy);
plot(xx1,yv)
ys=trapz(xx1,yv);
plot(xx2,ys)
yvd=derend(yv);
yvd=detrend(yv);
plot(xx1,yvd)
pv=polyfit(xx1,yv,1)
pv=polyfit(xx1',yv,1)
ppv=polyval(pv,xx1');
plot(xx1,yv-ppv)
hold on
plot(xx,yy)
plot(xx,yy,'r')
what
Fs = 1/0.1;  % Sampling Frequency
N   = 3;    % Order
Fc1 = 2;    % First Cutoff Frequency
Fc2 = 9;  % Second Cutoff Frequency
% Construct an FDESIGN object and call its BUTTER method.
h  = fdesign.bandpass('N,Fc1,Fc2', N, Fc1, Fc2, Fs);
Hd = butter(h);
Fs = 1/0.1;  % Sampling Frequency
N   = 4;    % Order
Fc1 = 2;    % First Cutoff Frequency
Fc2 = 8;  % Second Cutoff Frequency
% Construct an FDESIGN object and call its BUTTER method.
h  = fdesign.bandpass('N,Fc1,Fc2', N, Fc1, Fc2, Fs);
Hd = butter(h);
Fs = 1/0.1;  % Sampling Frequency
N   = 5;    % Order
Fc1 = 2;    % First Cutoff Frequency
Fc2 = 9;  % Second Cutoff Frequency
% Construct an FDESIGN object and call its BUTTER method.
h  = fdesign.bandpass('N,Fc1,Fc2', N, Fc1, Fc2, Fs);
Hd = butter(h);
Fs = 1/0.1;  % Sampling Frequency
N   = 3;    % Order
Fc1 = 2;    % First Cutoff Frequency
Fc2 = 5;  % Second Cutoff Frequency
% Construct an FDESIGN object and call its BUTTER method.
h  = fdesign.bandpass('N,Fc1,Fc2', N, Fc1, Fc2, Fs);
Hd = butter(h);
Fs = 1/0.1;  % Sampling Frequency
N   = 4;    % Order
Fc1 = 2;    % First Cutoff Frequency
Fc2 = 5;  % Second Cutoff Frequency
% Construct an FDESIGN object and call its BUTTER method.
h  = fdesign.bandpass('N,Fc1,Fc2', N, Fc1, Fc2, Fs);
Hd = butter(h);
yf=fiter(Hd,yy);
yf=filter(Hd,yy);
plot(xx,yf)
Fs = 1/0.1;  % Sampling Frequency
N   = 4;    % Order
Fc1 = 0;    % First Cutoff Frequency
Fc2 = 5;  % Second Cutoff Frequency
% Construct an FDESIGN object and call its BUTTER method.
h  = fdesign.bandpass('N,Fc1,Fc2', N, Fc1, Fc2, Fs);
Hd = butter(h);
Fs = 1/0.1;  % Sampling Frequency
N   = 4;    % Order
Fc1 = 0.1;    % First Cutoff Frequency
Fc2 = 5;  % Second Cutoff Frequency
% Construct an FDESIGN object and call its BUTTER method.
h  = fdesign.bandpass('N,Fc1,Fc2', N, Fc1, Fc2, Fs);
Hd = butter(h);
figure
yf=filter(Hd,yy);
plot(xx,yf)
yfv=trapz(xx,yf);
hold on
plot(xx1,yfv,'r')
yfvd=drend(yfv);
yfvd=derend(yfv);
yfvd=detrend(yfv);
plot(xx1,yfvd,'y);
plot(xx1,yfvd,'y');
yfs=trapz(xx1,yfvd);
figure
plot(xx2,yfs)
yfs=trapz(xx1,yfv);
hold
hold on
plot(xx2,yfs)
yfs=trapz(xx1,trapz(xx,ff));
yfs=trapz(xx1,trapz(xx,yy));
plot(xx2,yfs,'k')
title('sin  -> trapz trapz')
title('sin -> a ->v ->s trapz ,')
plot (yfs, 'DisplayName', 'yfs', 'YDataSource', 'yfs'); figure(gcf)
cc
what
load run16t
y=run16t_data(:,2);
yind=y>1.6e-4;
sy=sum(yind);
sy
plot(yind)
yind=y>0.00016;
sy=sum(yind)
plot(yind,'o')
3837/8192
plot(y)
x=run16t_data(:,1);
plot(x,y)
grid
yind=abs(y)>0.00016;\
yind=abs(y)>0.00016;
sum(yind)
plot(yind)
8192-7599
yind=abs(y)<0.00016;
yi=y(yind);
plot(yi)
%判断这些数据是否都是大于分辨率的
%分辨率的判断是要 用绝对值的，因为是幅值而不是代数值，加速度是可以有正有负的
title('resolutin  yind  yind=abs(y)>0.00016; and there is  yind=abs(y)<0.00016;')
cc
%-- 08-3-17 下午8:25 --%
what
%-- 08-3-17 下午9:14 --%
cc
up
u
up
cc
doc balance
%-- 08-3-18 下午2:40 --%
mex -setup
mbuild -setup
!cd e:\MATLAB7\work
!mccsavepath e:\MATLAB7\work
!mcc savepath e:\MATLAB7\work
test1(3)
mcc -t -L cpp -h test1.m
mcc -?
mcc -W lib -T test1.m
mcc -W  -T test1.m
mcc  -T test1.m
mcc  -m test1.m
doc mcc
mcc  -l test1.m
mcc -d target -m test.m
mcc -m -d target  test.m
mcc -d target -m example
mcc -W lib:libtriangle -T link:exe sierpinski.m triangle.c
mex -setup
mbuild -setup
mcc -W lib:libtriangle -T link:exe sierpinski.m triangle.c
test2
mcc -d target -m test2.m
cc
load count.dat
[n,p]=size(count)
plot(count)
t=1:n;
set(0,'defaultaxeslinestyleorder','-|--|-.')
set(0,'defaultaxescolororder',[0 0 0])
plot(t,count),
legend('Location 1','Location 2','Location 3',2)
xlabel('Time hour')
ylabel('Traffic count')
grid on
%-- 08-3-19  下午12:37 --%
doc inputdlg
doc tfdata
simulink
simout
doc cov
A = [-1 1 2 ; -2 3 1 ; 4 0 3]
cov(A)
diag(ans)
ans'
doc md
%-- 08-3-24 下午5:43 --%
help poly
ROOTS(POLY(1:20))
A =
1    2    3
4    5    6
7    8    0
A =...
1    2    3
4    5    6
7    8    0
A =...
[1    2    3
4    5    6
7    8    0]
p = poly(A)
roots(p)
eig(A)
z = eig(A)
c = zeros(n+1,1)
n=3
c = zeros(n+1,1)
c(1) = 1
for j = 1:n
c(2:j+1) = c(2:j+1)-z(j)*c(1:j)
end
%-- 08-3-26 下午5:54 --%
helpwin
doc quad
dblquad(@(x,y)sqrt(max(1-(x.^2+y.^2),0)), -1, 'x', -1, 1)
cc
demo simulink
plot (THETAdot, 'DisplayName', 'THETAdot', 'YDataSource', 'THETAdot'); figure(gcf)
plot (ForceF, 'DisplayName', 'ForceF', 'YDataSource', 'ForceF'); figure(gcf)
Y
clc
plot(Y)
plot (ForceF, 'DisplayName', 'ForceF', 'YDataSource', 'ForceF'); figure(gcf)
x0
plot (x, 'DisplayName','x', 'YDataSource', 'x'); figure(gcf)
plot (tt, 'DisplayName', 'tt', 'YDataSource', 'tt'); figure(gcf)
plot (t, 'DisplayName', 't', 'YDataSource', 't'); figure(gcf)
kr
kf
whos
h
clc
whos
plot(h)
cr,cr,kr,Iyy
Mb
Iyy
plot(t,h)
plot(tt,h)
L