[已答复] 解二元方程组,我编了一个程序无解,求大佬看一下问题出在哪了

[复制链接]
zgl123456654321 发表于 6 天前
想求出x、y、z但是运行显示无解
以下是程序代码:
clc; clear all;close all
E=116000;
h=25700;
v=0.3;
a=14.7*(10^-6);
k=0.2;
T=600;
K=450;
x1=175;%X方向应力
y1=193;%Y方向应力
z1=212;%Z方向应力
xz1=102;%XZ方向应力

axx=0;%第一步
ayy=0;%第一步
azz=0;%第一步
axz=0;%第一步

syms x z

Sxx=x1-1/3*(x1+y1+z1);%第一步
Syy=y1-1/3*(x1+y1+z1);%第一步
Szz=z1-1/3*(x1+y1+z1);%第一步
Sxz=xz1;%第一步
nxx=(Sxx-axx)/(2^0.5*K);%第一步
nyy=(Syy-ayy)/(2^0.5*K);%第一步
nzz=(Szz-azz)/(2^0.5*K);%第一步
nxz=(Sxz-axz)/(2^0.5*K);%第一步
eq1=1/E*(x-v*((1/2*(x+z))+z))+a*T+1/h*(x*nxx+(1/2*(x+z))*nyy+z*nzz+2*xz1*nxz)*nxx-(0.1587*((1/E*(x1-v*((1/2*(x+z))+z1)))+a*T+1/h*((x1*nxx+(1/2*(x+z))*nyy+z1*nzz+2*xz1*nxz)*nxx)));%第一个式子
eq2=1/E*((1/2*(x+z))-v*(x+z))+a*T+1/h*(x*nxx+(1/2*(x+z))*nyy+z*nzz+2*xz1*nxz)*nyy;%第二个式子
[x,z]=solve('eq1=0,eq2=0','x','z')
x=double(x);
z=double(z);
y=1/2*(x+z);


5 条回复


20141303 发表于 6 天前
仅供参考
  1. clc; clear all;close all
  2. E=116000;
  3. h=25700;
  4. v=0.3;
  5. a=14.7*(10^-6);
  6. k=0.2;
  7. T=600;
  8. K=450;
  9. x1=175;%X方向应力
  10. y1=193;%Y方向应力
  11. z1=212;%Z方向应力
  12. xz1=102;%XZ方向应力

  13. axx=0;%第一步
  14. ayy=0;%第一步
  15. azz=0;%第一步
  16. axz=0;%第一步

  17. syms x z

  18. Sxx=x1-1/3*(x1+y1+z1);%第一步
  19. Syy=y1-1/3*(x1+y1+z1);%第一步
  20. Szz=z1-1/3*(x1+y1+z1);%第一步
  21. Sxz=xz1;%第一步
  22. nxx=(Sxx-axx)/(2^0.5*K);%第一步
  23. nyy=(Syy-ayy)/(2^0.5*K);%第一步
  24. nzz=(Szz-azz)/(2^0.5*K);%第一步
  25. nxz=(Sxz-axz)/(2^0.5*K);%第一步
  26. eq1=1/E*(x-v*((1/2*(x+z))+z))+a*T+1/h*(x*nxx+(1/2*(x+z))*nyy+z*nzz+2*xz1*nxz)*nxx-(0.1587*((1/E*(x1-v*((1/2*(x+z))+z1)))+a*T+1/h*((x1*nxx+(1/2*(x+z))*nyy+z1*nzz+2*xz1*nxz)*nxx)));%第一个式子
  27. eq2=1/E*((1/2*(x+z))-v*(x+z))+a*T+1/h*(x*nxx+(1/2*(x+z))*nyy+z*nzz+2*xz1*nxz)*nyy;%第二个式子
  28. [x,z]=solve([eq1==0,eq2==0],[x,z])
  29. x=double(x);
  30. z=double(z);
  31. y=1/2*(x+z);
复制代码

zgl123456654321 发表于 6 天前

谢谢楼主,出来结果了,但是结果是两个数字相除,这个怎么改成具体的数值啊。我刚学matlab还不是很熟练

20141303 发表于 6 天前
仅供参考
  1. clc; clear all;close all
  2. E=116000;
  3. h=25700;
  4. v=0.3;
  5. a=14.7*(10^-6);
  6. k=0.2;
  7. T=600;
  8. K=450;
  9. x1=175;%X方向应力
  10. y1=193;%Y方向应力
  11. z1=212;%Z方向应力
  12. xz1=102;%XZ方向应力

  13. axx=0;%第一步
  14. ayy=0;%第一步
  15. azz=0;%第一步
  16. axz=0;%第一步

  17. syms x z

  18. Sxx=x1-1/3*(x1+y1+z1);%第一步
  19. Syy=y1-1/3*(x1+y1+z1);%第一步
  20. Szz=z1-1/3*(x1+y1+z1);%第一步
  21. Sxz=xz1;%第一步
  22. nxx=(Sxx-axx)/(2^0.5*K);%第一步
  23. nyy=(Syy-ayy)/(2^0.5*K);%第一步
  24. nzz=(Szz-azz)/(2^0.5*K);%第一步
  25. nxz=(Sxz-axz)/(2^0.5*K);%第一步
  26. eq1=1/E*(x-v*((1/2*(x+z))+z))+a*T+1/h*(x*nxx+(1/2*(x+z))*nyy+z*nzz+2*xz1*nxz)*nxx-(0.1587*((1/E*(x1-v*((1/2*(x+z))+z1)))+a*T+1/h*((x1*nxx+(1/2*(x+z))*nyy+z1*nzz+2*xz1*nxz)*nxx)));%第一个式子
  27. eq2=1/E*((1/2*(x+z))-v*(x+z))+a*T+1/h*(x*nxx+(1/2*(x+z))*nyy+z*nzz+2*xz1*nxz)*nyy;%第二个式子
  28. [x,z]=vpasolve([eq1==0,eq2==0],[x,z])
  29. x=double(x);
  30. z=double(z);
  31. y=1/2*(x+z);
复制代码

zgl123456654321 发表于 4 天前

感谢,又学到一个新知识:lol

zgl123456654321 发表于 4 天前

还想在请教一个问题。给定x,y,z的初始值,他们每次各自减去各自的固定值,条件j大于0的时候,用上一次的差值,再次减去那个固定值;当条件j小于等于0的时候,输出结果x,y,z。所编的代码如下:
clc; clear all;close all
%设定初始参数
x=2356;%初始应力数值
y=2768;
z=1952;
j=1;%初始迭代条件
cs=1;%初始迭代次数
while j>0     %设定迭代判断条件数值
syms x y z
%更新三个方向上的应力
x=x-532;
y=y-551;
z=z-570;
%更新迭代判定条件
j=1/(2^(0.5))*((x-y)^2+(y-x)^2+(z-x)^2)-606^2;
cs=cs+1;
end
您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
热门教程
站长推荐
快速回复 返回顶部 返回列表