solve,subs循环求解方程组

%本程序是求解原点O(0,0)与(X,Y)的连线L1与AB对应17条直线的交点的问题
%也就是要求出17个交点
clear,clc
%此处A和B分别对应18个点
A=[10 20 30 40 50 60 67 68 67.5 65 60 58 50 40 30 20 18 10];
B=[88 84 77 75 82 86 80 70 60 50 44 40 35 36 40.5 48 50 62];
%plot(A,B)
X=25.645;
Y=61.774;
%两个点:点1(A_i,B_i),
%点2(A_i_1,B_i_1)
syms A_i B_i A_i_1 B_i_1 x y X1 Y1                                                %此处的定义很重要
S=solve('(B_i_1-B_i)*x-(A_i_1-A_i)*y+((A_i_1-A_i)*B_i_1-(B_i_1-B_i)*A_i_1)=0',...     %方程1
        'Y1*x-X1*y=0','x','y');                                                                   %方程2
    S.x,S.y
    for i=1:17;
    m=S.x;
    n=S.y;
    A_ii=A(i);
    B_ii=B(i);
    A_ii_1=A(i+1);
    B_ii_1=B(i+1);
    m=subs(m,A_i,A_ii);                         %subs(m,a,b),表示用b将m中的a替换,m为关于a的符号表达式
    m=subs(m,B_i,B_ii);                         %如m=a+5+c,执行subs(m,a,b)后:m=b+5+c;
    m=subs(m,A_i_1,A_ii_1);                 %其中b也可以是数字,是符号函数向数值计算的一种转化.
    m=subs(m,B_i_1,B_ii_1);                
    n=subs(n,A_i,A_ii);
    n=subs(n,B_i,B_ii);
    n=subs(n,A_i_1,A_ii_1);
    n=subs(n,B_i_1,B_ii_1);
    m=subs(m,X1,X);
    m=subs(m,Y1,Y);
    n=subs(n,X1,X);
    n=subs(n,Y1,Y);
    data(1,i)=m;
    data(2,i)=n;
end
data
plot(data(1,:),data(2,:),'-*')
%上面的程序到此结束

 

%subs用法个人心得:
clear,clc
x=0:0.01:2*pi;
y=sin(x);
subplot(3,1,1);plot(x,y);
title('原图');
%--------------------
syms t y1;
y1='sin(t)';
subplot(3,1,2);ezplot(y1);
title('符号函数画图');
%--------------------
n=0:0.01:3*pi;
y2=subs(y1,t,n);      %表示将y1=sin(t)中的t替换为n,
subplot(3,1,3);plot(n,y2);           %也就是将符号变换为数值
title('替换后的绘图');

%函数文件

%求两条直线的交点,两条直线为
% L1: A1(m1,n1),B1(m2,n2)所确定的直线 A1*x+B1*y+C1=0;
% L2: A2(m3,n3),B2(m4,n4)所确定的直线 A2*x+B2*y+C2=0;
function [Pt]=solve_point(m1,n1,m2,n2,m3,n3,m4,n4)
A1=n2-n1; B1=m1-m2; C1=m2*n1-m1*n2;
A2=n4-n3; B2=m3-m4; C2=m4*n3-m3*n4;
A=[A1,B1;A2,B2];B=[-C1;-C2];
Pt=A/B;

 

 

 

 

 

%改进

tic
%本程序是求解原点O(0,0)与(X,Y)的连线L1与AB对应17条直线的交点的问题
%也就是要求出17个交点
clear,clc
%此处A和B分别对应18个点
A=[10 20 30 40 50 60 67 68 67.5 65 60 58 50 40 30 20 18 10];
B=[88 84 77 75 82 86 80 70 60 50 44 40 35 36 40.5 48 50 62];
plot(A,B,'-*')
X=25.645;
Y=61.774;
hold on;plot([0,40],[0,96.3525],'-r');
Pt=solve_point(0,0,X,Y,30,77,40,75)
text(30,84,'(31.8152,76.6370)');
for i=1:17;
    [Pt]=solve_point(0,0,X,Y,A(1,i),B(1,i),A(1,i+1),B(1,i+1));
    data(:,i)=Pt;
end
data'
toc