Matlab实验（一）

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)