Matlab实验(一)
1、下面给出一个迭代模型
写出求解该模型的M-函数。如果迭代初值为 ,那么请进行30000次迭代求出一组x和y向量,然后再所有的 坐标处点亮一个点(注意不要连线),最后绘制出所需的图形。(注:这样绘制的图形又称为Henon引力线图,它将迭代出来的随机点吸引到一起,最后得出貌似连贯的引力线图。)
function [x,y]=Henons(x0,y0,n)
%x0、y0:初始值
%n:迭代次数
%返回相应的x,y的坐标值
x(1)=x0;
y(1)=y0;
for i=2:n
x(i)=1+y(i-1)-1.4*x(i-1)*x(i-1);
y(i)=0.3*x(i-1);
end
plot(x,y,'*');
grid on;
title('Henon引力线图');
end
%调用
[x,y]=Henons(1,1,30000)
2.试用两种方法求解代数方程 。
(1) 二分法。若在某个区间(a,b)内,f(a)*f(b)<0 ,则该区间存在方程的根。取中点 ,则可以根据 f(x1)和 f(a)、f(b)的关系确定根的范围,用这样的方法可以将区间的长度减半。重复这样的过程,直至区间长度小于预先指定的ε,则可以认为得出的区间端点是方程的解。令 ,试用二分法求区间(-4,0)内方程的解。
(2) 梯度法。假设该方程解的某个初始猜测点为 xn,则梯度法可以得出下一个近似点 。若两个点足够近,即 ,其中ε为预指定的误差限,则认为 是方程的解,否则将 设置为初值继续搜索,直到得出方程的解。令 x0=-4,,试用梯度法求解上面的方程。
(1) 二分法
%%%%%%方法一:function版本
function answer=Binary_solve(f,x,a,b,precision)
%该算法由符号变量实现
%f:相应的方程 x:自变量
%a:左区间 b:右区间
%precision:准确度
if subs(f,x,a)*subs(f,x,b)>0
disp('条件不满足');
return
else
while(abs(b-a)>precision)
mid=(b+a)/2;
if (subs(f,x,a)*subs(f,x,mid))<0
b=mid;
elseif (subs(f,x,a)*subs(f,x,mid)>0)
a=mid;
else
answer=mid;%中等于0,提前退出循环
break;
end
end
answer=mid;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
syms f x;
f=x*x*sin(0.1*x+2)-3;
root=Binary_solve(f,x,-4,0,10^(-10))
root =
-1.7603
%%%%%%%%%%%方法二
%%%%%%%%%%%
f=@(x)(x.^2*sin(0.1.*x+2)-3)
a=-4;
b=0
root=[]
while abs(b-a)>1e-10
temp=(b+a)./2;
if f(b).*f(temp)<0
a=temp;
elseif f(a).*f(temp)<0
b=temp;
else
if f(b)==0
root=[root,b];
elseif f(a)==0
root=[root,a];
else
root=[root,temp];
end
break;
end
end
root=[root,temp]
(2) 梯度法
syms f x
f=x*x*sin(0.1*x+2)-3;
x0=-4;
prec=1e-12;
dif=diff(f);
x1=x0;
x2=x0;
while(1)
x1=x2;
x2=x1-subs(f,x,x1)/subs(dif,x,x1);
if abs(x2-x1)<prec
break;
end
end
answer=double(x2)
3.根据 ,求 的近似值。当n分别取100,1000,10000时,结果是多少?
(要求:分别用循环结构和向量运算(使用sum函数)来实现)
function [answer1,answer2]=pi2_div_6(n)
%answer1:sum实现pi……2/6
%answer2:for循环实现pi/6
i=1:n;
answer2=0;
disp('sum实现计时');
tic;
answer1=sqrt(6*sum(1./i./i));
toc;
disp('for循环实现计时');
tic;
for i=1:n
answer2=answer2+1/i/i;
end
answer2=sqrt(6*answer2);
toc;
end
4.根据 ,求
(1) y<3时的最大n值;
(2) 与(1)的n值对应的y值。
i=1;
y=0;
while(1)
y=y+1/(2*i-1);
if y>3
break;
end
i=i+1;%注意不要写到上面的end种
end
answer1=i-1;
answer2=y-1/(2*i-1);
end
5. 已知某迭代序列 ,并已知该序列当n足够大时将趋于某个固定的常数,试选择合适的n值,找出该序列的稳态值(达到精度要求10-14),并找出其精确的数学表示。
x0=1;
prec=1e-14;
i=2;
x(1)=x0;
x(i)=x(i-1)/2+3/2/x(i-1);
while(1)
x(i)=x(i-1)/2+3/2/x(i-1);
if(x(i)-x(i-1)<prec)
break;
end
%实际上需要几个连续的相邻项之差都要小于prec才能达到稳态
i=i+1;
end
x=x(i);
6.
function y=twof_solve( Cauchy,left,right,temp)
%Cauchy:阻尼比(行向量) left左边区间 right:右边区间
%temp:指定间隔默认0.1,越小曲线越精确
%y:返回给定阻尼比的响应值
syms t f;%Bera由于是系统变量,所以命名实际函数
k=left:temp:right;%待测区间
figure(1);
hold on;
grid on;
y=zeros(length(Cauchy),length(k));
legends='';
for i=1:length(Cauchy)
%每一个Cauchy 都对于一个不同的函数,所以需要在循环内写函数
Beta=sqrt(1-Cauchy(i)*Cauchy(i));
Theta=atan(sqrt(1-Cauchy(i)*Cauchy(i))/Cauchy(i));
f=1-1/Beta*exp(-Cauchy(i)*t)*sin(Beta*t+Theta);
y(i,:)=subs(f,t,k);
plot(k,y(i,:),'-');
legends=[legends;['Cauchy=',num2str(Cauchy(i))]];
end
legend(legends);%legend不能写在循环内,否者会被覆盖
title('不同ε的响应曲线图');
xlabel('t');
ylabel('y');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%调用
twof_solve([0.2 0.4 0.6 0.8],0,18,0.01)
图像:
————————————————————————————————————————————————————————————————————————————————
说明:以上习题来自《控制系统计算机辅助设计 —— MATLAB语言与应用( 第3版)》,薛定宇老师的著作,解答仅供参考,如有错误,请指正。
以下是第二章习题的一些解答,仅供参考。
%2-5
syms x s;
f=x*sin(x)/((x^2+2)^0.5*(x+5));
pretty(f)
x=(s-1)/(s+1);
pretty(f)
pretty(subs(f,x))
%2-7
%数值解
clear
tic;
n=63;
s=ones(n,1);
for i=2:n
%s(i)=2^(i-1);
s(i)=2*s(i-1);
end
disp('数值解=')
disp(sum(s))
toc;
tic;
n=63;
s=ones(n,1);
sym(s);
for i=2:n
%s(i)=2^(i-1);
s(i)=2*s(i-1);
end
disp('符号解=')
disp(sum(s))
toc;
%2-9(寻找质数)
%使用封装好的库函数;
clear;
n=1000;
pri=primes(n);
%手撸代码
prime_num=[]
for i=2:n
k=1;
i
for j=2:fix(sqrt(i))
k=k+1
if(mod(i,j)==0)
k=k-1
break;
end
end
if k>=fix(sqrt(i))
prime_num=[prime_num,i];
i
end
end
pri==prime_num
%2-10
syms n s
s=0;
n=1;
E=1;
while(abs(E)>=1e-12)
s=s+1+2/(n*n);
% if(fix(s)==s)
% n=n+1;
% continue;
% end
% E=s-fix(s)
E=1/(n*n)
n=n+1;
end
disp(s)
disp(n)
%2-14
x=input('请输入x的值:');
h=input('请输入h的值:');
D=input('请输入D的值:');
if x>D
f=h
elseif abs(x)<=D
f=h/D*x
elseif x<-D
f=-h
end
%2-22
syms t;
t1=-1:0.1:0;
t2=0:0.1:1;
y1=sin(ones(length(t1)-1)/t1(1:end-1));
y2=sin(ones(length(t2)-1)/t2(2:end));
t0=0;
y0=double(limit(sin(1/t),0))
t=[t1(1:end-1),t0,t2(2:end)];
y=[y1',y0,y2']
plot(t,y)
%2-24
theta = 0:0.01:2*pi;
figure(1);
p1=1.0013.*theta.^2;
p2=cos(7.*theta./2);
p3=sin(theta)./theta;
p4=1-cos(7.*theta).^3;
subplot(421);
polar(theta,p1);
subplot(422);
polar(theta,p2);
subplot(423);
polar(theta,p3);
subplot(424);
polar(theta,p4);
%2-26
[x,y] = meshgrid([-5:0.1:5],[-5:0.1:5]);
z1=x.*y;
z2=sin(x.*y);
z3=sin(x.*x+y.*y);
z4=-x.*y.*exp(-2.*(x.^2+y.^2));
figure(1);
subplot(421);
surf(x,y,z1)
subplot(422);
surfc(x,y,z1)
subplot(423);
surfl(x,y,z1)
subplot(424);
waterfall(x,y,z1)
figure(2);
subplot(421);
surf(x,y,z2)
subplot(422);
surfc(x,y,z2)
subplot(423);
surfl(x,y,z2)
subplot(424);
waterfall(x,y,z2)
figure(3);
subplot(421);
surf(x,y,z3)
subplot(422);
surfc(x,y,z3)
subplot(423);
surfl(x,y,z3)
subplot(424);
waterfall(x,y,z3)
figure(4);
subplot(421);
surf(x,y,z4)
subplot(422);
surfc(x,y,z4)
subplot(423);
surfl(x,y,z4)
subplot(424);
waterfall(x,y,z4)