查看: 173|回复: 3|关注: 0

[已答复] 求帮忙改程序,将下列程序由求Za Zb Zc Ze解析解改为求数值解。

[复制链接]

新手

10 麦片

财富积分


050


1

主题

3

帖子

0

最佳答案
发表于 2020-7-4 11:20:12 | 显示全部楼层 |阅读模式
%% 以下是产生空间四面体ABCE空间坐标的程序
clc;close all;
global f L1 L2 x_a x_b x_c x_e y_a y_b y_c y_e;
xe=1;                                     %输入E1点x坐标
ye=1;                                       %输入E1点y坐标
ze=0;                                       %输入E1点z坐标
r=10;                                        %输入球面半径
L1=sqrt(3)*r;
f=1;                                         %输入焦距
n=[1,0,0];                                   %三角形面法线n
n0=1;                                        %三角形顶点参数
xed=xe+n0*(n(1)/sqrt(n(1)^2+n(2)^2+n(3)^2)*r);     %产生四面体顶点坐标
yed=ye+n0*(n(2)/sqrt(n(1)^2+n(2)^2+n(3)^2)*r);
zed=ze+n0*(n(3)/sqrt(n(1)^2+n(2)^2+n(3)^2)*r);     %产生四面体顶点坐标
L2=sqrt((xed-xe)^2+(yed-ye)^2+(zed-ze)^2+r^2);
fenzhibianliang=0;                          %该分支变量用于选择显示哪个三角形,
syms xA yA zA                               %以E点x值为A点x值,解出符合条件的A点两个,取其中一点。
eq1 = (xA-xe)^2+(yA-ye)^2+(zA-ze)^2==r^2;
eq2 = n(1) *(xA-xe)+ n(2) *(yA-ye)+ n(3) *(zA-ze)==0;%以上的是生成圆环
eq3 = xA==xe;
[xA, yA, zA] = solve(eq1, eq2, eq3, xA, yA, zA);
m=[(xA(2)-xe)/r, (yA(2)-ye)/r, (zA(2)-ze)/r];                    %得出第一个切割面向量m

syms xB yB zB
eq1 = (xB-xe)^2+(yB-ye)^2+(zB-ze)^2==r^2;
eq2 = n(1) *(xB-xe)+ n(2) *(yB-ye)+ n(3) *(zB-ze)==0;%以上的是生成圆环
eq3 = m(1)*(xB-xe)/r + m(2)*(yB-ye)/r + m(3)*(zB-ze)/r==-0.5;
[xB, yB, zB] = solve(eq1, eq2, eq3, xB, yB, zB);
m1=[(xB(1)-xe)/r, (yB(1)-ye)/r, (zB(1)-ze)/r];                  %得出第二个切割面向量m1
m2=[(xB(2)-xe)/r, (yB(2)-ye)/r, (zB(2)-ze)/r];                  %得出第三个切割面向量m2

A=[(xA(2)),(yA(2)),(zA(2))]                                          %确定并输出ABCE空间坐标
B=[(xB(1)),(yB(1)),(zB(1))]
C=[(xB(2)),(yB(2)),(zB(2))]
E=[(xed),(yed),(zed)]
%% 以下是用来计算平面投影坐标的

X_a=A(1);
Y_a=A(2);
Z_a=A(3);
X_b=B(1);
Y_b=B(2);
Z_b=B(3);
X_c=C(1);
Y_c=C(2);
Z_c=C(3);
X_e=E(1);
Y_e=E(2);
Z_e=E(3);
x_a=f*X_a/(f-Z_a);
y_a=f*Y_a/(f-Z_a);
x_b=f*X_b/(f-Z_b);
y_b=f*Y_b/(f-Z_b);
x_c=f*X_c/(f-Z_c);
y_c=f*Y_c/(f-Z_c);
x_e=f*X_e/(f-Z_e);
y_e=f*Y_e/(f-Z_e);
%% 以下是反推空间坐标的
syms Za Zb Zc Ze
eq1 = (Za*x_a/f-Zb*x_b/f)^2+(Za*y_a/f-Zb*y_b/f)^2+(Za-Zb)^2==L1^2;
eq2 = (Za*x_a/f-Zc*x_c/f)^2+(Za*y_a/f-Zc*y_c/f)^2+(Za-Zc)^2==L1^2;
eq3 = (Zb*x_b/f-Zc*x_c/f)^2+(Zb*y_b/f-Zc*y_c/f)^2+(Zb-Zc)^2==L1^2;
eq4 = (Ze*x_e/f-Za*x_a/f)^2+(Ze*y_e/f-Za*y_a/f)^2+(Ze-Za)^2==L2^2;
eq5 = (Ze*x_e/f-Zb*x_b/f)^2+(Ze*y_e/f-Zb*y_b/f)^2+(Ze-Zb)^2==L2^2;
eq6 = (Ze*x_e/f-Zc*x_c/f)^2+(Ze*y_e/f-Zc*y_c/f)^2+(Ze-Zc)^2==L2^2;
[Za,Zb,Zc,Ze] = solve(eq1, eq2, eq3, eq1, eq4, eq5, eq6, Za,Zb,Zc,Ze);

回复主题 已获打赏: 0 积分

举报

新手

10 麦片

财富积分


050


1

主题

3

帖子

0

最佳答案
 楼主| 发表于 2020-7-4 11:21:03 | 显示全部楼层
请大佬指点
回复此楼 已获打赏: 0 积分

举报

论坛优秀回答者

7

主题

1508

帖子

321

最佳答案
  • 关注者: 77
发表于 2020-7-4 11:47:19 | 显示全部楼层
最后一句改成
  1. [Za,Zb,Zc,Ze] = vpasolve([eq1, eq2, eq3, eq1, eq4, eq5, eq6], [Za,Zb,Zc,Ze]);
复制代码

试试
回复此楼 已获打赏: 0 积分

举报

新手

10 麦片

财富积分


050


1

主题

3

帖子

0

最佳答案
 楼主| 发表于 2020-7-4 16:18:03 | 显示全部楼层
20141303 发表于 2020-7-4 11:47
最后一句改成

试试

并不是将解析解改为数值解,而是从一开始解最后的方程组时就是求出数值解。我现在的问题是只要代码开头的预设条件变的比较不规整,最后一个方程组解析解就求不出来。我想要的是利用数值分析的方法求解最后一个方程组的代码。谢谢大佬解惑。
回复此楼 已获打赏: 0 积分

举报

您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

站长推荐上一条 /4 下一条

快速回复 返回顶部 返回列表