哈工大数学实验报告(完整版).pdf

实验一 Matlab 的使用 1. 输入 A,B A = -3 5 0 8; 1 -8 2 -1;0 -5 9 3; -7 0 -4 5 A = -3 5 0 8 1 -8 2 -1 0 -5 9 3 -7 0 -4 5 B = 0 2 -1 6 B = 0 2 -1 6 (a) 方法 AB ans = -0.6386 -0.4210 -0.3529 0.0237 (b) 方法 inv(A)*B ans = -0.6386 -0.4210 -0.3529 0.0237 2. 分块矩阵 E = eye(3); R = rand(3,2); O = zeros(2,3); S = diag(1 2); A = E R;O S A = 1.0000 0 0 0.8147 0.9134 0 1.0000 0 0.9058 0.6324 0 0 1.0000 0.1270 0.0975 0 0 0 1.0000 0 0 0 0 0 2.0000 A2 = E R+R*S;O S2; err = A2 - A2 norm(err) ans = 0 3. 多项式 r = 2 -3 1+2*i 1-2*i 0 -6; p = poly(r) p = 1 5 -9 -1 72 -180 0 计算 0.8-1.2 处的值 y1 = polyval(p,0.8) y1 = -100.2179 y1 = polyval(p,-1.2) y1 = 293.2900 4. 特征方程和友元矩阵 p = 1 0 -6 3 -8; A = compan(p) A = 0 6 -3 8 1 0 0 0 0 1 0 0 0 0 1 0 求 A 的特征值 lamda = eig(A) lamda = -2.8374 2.4692 0.1841 + 1.0526i 0.1841 - 1.0526i roots(p)结果 roots(p) ans = -2.8374 2.4692 0.1841 + 1.0526i 0.1841 - 1.0526i 结论一维向量 p 友元矩阵的特征值即为其降幂多项式的根 5. 做函数的图形 x = linspace(0,1); X = x.2;x.3;x.4;x.5; plot(x, X(1,:), r, x, X(2,:), g, x, X(3,:), b, x, X(4,:), k) legend(y = x2,y = x3,y = x4,y = x5) 用 subplot 画四个函数图形的指令 subplot(2,2,1); plot(x,X(1,:); title(y = x2); subplot(2,2,2); plot(x,X(2,:); title(y = x3); subplot(2,2,3); plot(x,X(3,:); title(y = x4); subplot(2,2,4); plot(x,X(4,:); title(y = x5); 6. fplot 的使用 fplot(funfplot,-0.1 0.1, 2e-4) 7. 用作图法求 的根的近似值 x = linspace(-6,4); y = 4*sin(x)-x-2; plot(x,y) y0,n = find(abs(y) plot(x, y, x(n), y(n), ro) 对 x,y 轴加标志并加网格线 xlabel(x axis) ylabel(y axis) grid on 8. 作曲面 的三维图形 X,Y = meshgrid(-5:0.1:5); Z = X.2-Y.2; meshc(Z) 或使用 surfc(Z) 9. 建立 m 函数 (a) 自然数阶乘 function y = Factorial(x) if round(x) = x | x 0 error(x 必须是非负整数); end if x = 0 y = 1; else y = prod(1:x); end (b) n 中取 m 的组合 function y = Combination(n,m) if round(m) = m | round(n) = n error(m,n 必须为整数); elseif n 1 | m 0 | n taylor0246 input your function (only x variable): 1/720*x6 - 2*x2 +exp(x) input the point to use taylor: 1 实验二 平面线性映射的迭代 1. 求特征值和特征向量 A = 3/4 1/2 1/4;1/8 1/4 1/2;1/8 1/4 1/4 A = 0.7500 0.5000 0.2500 0.1250 0.2500 0.5000 0.1250 0.2500 0.2500 V,D = eig(A) V = -0.9094 -0.8069 0.3437 -0.3248 0.5116 -0.8133 -0.2598 0.2953 0.4695 D = 1.0000 0 0 0 0.3415 0 0 0 -0.0915 2. 计算公式 n = 100 时计算结果为 P = 0 1/2 1/2; P100 = A100*P P100 = 0.6087 0.2174 0.1739 3. 与 A 的特征向量比较 P100/-0.6693-V(:,1) ans = 1.0e-004 * -0.5342 -0.1908 -0.1526 发现经过多次迭代后绝对值小于 1 的特征值对结果影响会越来越弱而绝对值等于一的 特征值对结果的影响不会变小 4. 求 A 特征值和特征向量 V,D = eig(A) V = 0.3827 -0.9239 -0.9239 -0.3827 D = -1.4142 0 0 1.4142 5. 点列轨迹 function IterationDraw A = 1 1;1 -1; V2 = 2;0; V3 = 2;2; V4 = 0;2; for i = 1:15 V2(:,i+1) = A*V2(:,i); V3(:,i+1) = A*V3(:,i); V4(:,i+1) = A*V4(:,i); end subplot(2,2,1); plot(V2(1,:),V2(2,:),go,V3(1,:),V3(2,:), b*,V4(1,:),V4(2,:),r+) subplot(2,2,2); plot(V2(1,:),V2(2,:),go) subplot(2,2,3); plot(V3(1,:),V3(2,:),b*) subplot(2,2,4); plot(V4(1,:),V4(2,:),r+) 6. 迭代映射 A = 1 2;3 -1 A = 1 2 3 -1 Planer(A,1 1,6) 7. Planar线性迭代 function Planar1(A, x0, maxit) m,n = size(A); x0 = x0(:); if m = 2 | n = 2 error(A must be a 2-dimention square matrix); elseif length(x0) = 2 error(A and x0 must agree); end X = zeros(length(x0), maxit+1); X(:,1) = x0; for i = 1:maxit X(:,i+1) = A*X(:,i); end plot(X(1,:), X(2,:), ro); Planar1(1/5 99/100;1 0,1 1,100) Planar1(1/5 99/100;1 0,1 0,100) 8. a的含义a为模最大的特征值 9. 经过多次迭代后模最大的特征值对结果的影响掩盖了其他特征值对结果的影响 10. Planar线性迭代 Planar1(0 1;100/99 -20/99,1 1,100) Planar1(0 1;100/99 -20/99,1 0,100) V,D = eig(0 1;100/99 -20/99) V = 0.7399 -0.6690 0.6727 0.7433 D = 0.9091 0 0 -1.1111 模最大的特征值 -1.1111属于这一特征值的方向-0.6690, 0.7433 11. 动画演示程序 function AnimDemo(a, b, alpha, beta, shape) e1 = cos(alpha);sin(alpha); e2 = cos(beta);sin(beta); A = a*e1, b*e2 * inv(e1, e2); if strcmp(shape,square) = 1 X0 = 1 -1 -1 1 1; 1 1 -1 -1 1; elseif strcmp(shape,circle) = 1 theta = linspace(0,2*pi); X0 = cos(theta);sin(theta); else error(Please decide shape); end X = A*X0; plot(X0(1,:), X0(2,:), r, X(1,:), X(2,:), b-) AnimDemo(1.2,0.8,pi/4,3*pi/4,square) 12. 变换圆形 theta = linspace(0,2*pi); x = cos(theta); y = sin(theta); Data = x;y; C = 1 0.2;0.2 1; Data1 = C*Data; plot(Data(1,:),Data(2,:),r,Data1(1,:),Data1(2,:),b) C的特征值和特征向量 V,D = eig(C) V = -0.7071 0.7071 0.7071 0.7071 D = 0.8000 0 0 1.2000 设计坐标变化矩阵,压缩拉伸方案 先确定拉伸、压缩的方向则方向向量 确定压缩拉伸比例,通过 13. 使用fill绘图函数 function fillTest x = 1 -1 -1 1 1; y = 1 1 -1 -1 1; k = 10; axis(-k k -k k); hold on fill(x,y,y); a1 = 1 1; a2 = -1 1; k1 = 1.2; k2 = 0.8; A = k1*a1 k2*a2*inv(a1 a2); Z = x;y; s = gbrygbrygb; for n = 1:10 Z = A*Z; pause(1) fill(Z(1,:),Z(2,:),s(n); end hold off pause(2); syms o R = cos(o) -sin(o);sin(o) cos(o); fill(Z(1,:),Z(2,:),y); Z0 = Z; hold on for theta = 0:0.1*pi:pi o = theta; A = eval(R); Z = A*Z0; pause(1); fill(Z(1,:),Z(2,:),y); end hold off 实验三 常微分方程数值解 2. 单摆问题 t 的维度与 ode23 函数计算时选择的自适应步长有关步长越小t 维度越高 t 是 n 维列向量x 是 n*2 的矩阵第一列是的值第二列是的值 去掉第 6 行分号后xdot 的结果出现了 n 次n 是 t 的纬度 说明 ode23 这个函数每 步计算都要调用一次被求函数 t 在第七行语句中赋值 3. 改进欧拉公式 function x,y = EulerExODE(fun, xEnd, ini, h) ini = ini(:); x0 = ini(1); y0 = ini(2:length(ini); x = (x0:h:xEnd); y = ones(length(x),length(ini)-1); for i = 1:length(y0) y(:,i) = y(:,i)*y0(i); end for j=2:length(x) y(j,:) = y(j-1,:) + h*feval(fun, x(j-1), y(j-1,:); y(j,:) = y(j-1,:) + 0.5*h*( feval(fun, x(j-1), y(j-1,:) + feval(fun, x(j), y(j,:) ); end 或 x,y = EulerExODE(ode1, 2, 0 0); x1,y1 = ode23(ode1,0 2,0); plot(x,y,ro,x1,y1) x,y = EulerExODE(ode1, 2, 0 1); x1,y1 = ode23(ode1,0 2,1); plot(x,y,ro,x1,y1) x,y = EulerExODE(ode2, 2, 0 1 0); x1,y1 = ode23(ode2, 0 2, 1 0); plot(x,y,r,x1,y1,b*) 4. 小船过河问题 function dy, isterminal, direction = cross(t, y, flag) v1 = 1; v2 = 2; y0 = 0; 100; h = y0-y; nh = norm(h); if nargin CompeteTest(1,1,100,100,0.5,2,10,10) CompeteTest(1,1,100,100,0.5,2,10,5) CompeteTest(1,1,100,100,0.5,2,5,10) CompeteTest(2,1,100,100,0.5,2,10,10) CompeteTest(1,2,100,100,0.5,2,10,10) CompeteTest(1,1,100,100,1.5,0.7,10,10) CompeteTest(1,1,100,100,0.8,0.7,10,10) CompeteTest(1,1,100,100,1.5,1.7,10,10) 分析改变对整体趋势影响不大 影响 x,y 改变的速度并且越大系统到达稳定状态越快 影响系统的稳态若它们分布于 1 的两侧则 s 大于 1 的种群会淘汰另 一个种群若它们都小于 1稳态会是一个平衡共生的局面若它们都大于 1则 s 比较大的种群将淘汰另一个种群。 6. 导弹追踪飞机问题 function S = A(t) if t 6 S = 8*t;0; else S = 8*(12-t);1; end function zs,isterminal,direction = B(t,z,flag) global w X = A(t); h = X - z; nh = norm(h); if nargin 3 | isempty(flag) zs = (w/nh)*h; else switch(flag) case events zs = nh - 1e-3; isterminal = 1; direction = 0; otherwise error(unknown flag:, flag); end end function cross_1(cx,cy) v = 2; kx = cx,cx,cx,cx-v,cx+v; ky = cy,cy+2.5*v,cy+1.5*v,cy+1.5*v,cy+1.5*v; plot(kx,ky); plot(cx,cy,o); %main1.m global w yo = 60;70; w = 10; options = odeset(RelTol,1e-5,Events,on); t,y = ode23(B,0 20, yo, options); j = ; for h = 1:length(t) w = A(t(h); j = j;w; end xmin = min(min(y(:,1), min(j(:,1); xmax = max(max(y(:,1), max(j(:,1); ymin = min(min(y(:,2), min(j(:,2); ymax = max(max(y(:,2), max(j(:,2); clf; hold on; axis(xmin xmax ymin-1 ymax); for h = 1:length(t)-1 plot(y(h,1),y(h+1,1),y(h,2),y(h+1,2),-,. color,red,EraseMode,none); plot(j(h,1),j(h+1,1),j(h,2),j(h+1,2),:,. color,black,EraseMode,none); drawnow; pause(0.3); end 7. p = max(size(y); 8. cross_1(y(p,1),y(p,2); 9. hold off; 实验四 Matlab 程序编写 一、 options = odeset()的设定形式 1. 高阶微分方程转化为一阶微分方程组 等价的一阶微分方程组 2. 求解 Vandor Pol 方程 等价的一阶微分方程组 待求解的 B 的追踪路线程序 function ydot = vdpol(t,y,mu) if nargin 3 mu = 2; end ydot = y(2); mu*(1-y(1)2)*y(2)-y(1); 终止ode45调用的事件函数 function value, isterminal, direction = vdpolevents(t,y,mu) value(1) = abs(y(2) - 2; isterminal(1) = 0; direction(1) = 0; 求解方程的主程序 %mainvdpol.