필터 지우기
필터 지우기

Getting imaginary values as solutions to equations using solve function

조회 수: 1 (최근 30일)
AD
AD 2023년 8월 28일
댓글: Dyuman Joshi 2023년 8월 28일
I am trying to solve these three equations in matlab. I am getting the answer as a function of z and imaginary values after substituting the values of z and a warning. Where am I going wrong..can I get only real valued solutions? Also, I want to plot sol1 as a function of z. Can someone please help me with this?
P_l=50;
v=0.1;
k=15;
Tm=1375;
T0=300;
alpha=3.75*10^(-6);
E= 190*10^9;
h= 100*10^9;
G= 150*10^9;
nu=0.3;
psi= 1- exp(-(0.45*h)/(2*G));
syms x z x_prime z_prime t dT_dx dT_dz sigma_xx sigma_yy sigma_zz;
xm=x-x_prime;
zp=z+z_prime;
zm=z-z_prime;
% Define the terms
T=(P_l*exp(-v*((sqrt((x-v*t)^2 + (z)^2) + (x-v*t)))/(2*alpha)))/(4*3.14*k*sqrt((x-v*t)^2 + (z)^2)) +T0;
dT_dx=diff(T,x);
dT_dx_prime=subs(dT_dx,[x,z],[x_prime,z_prime]);
dT_dz=diff(T,z);
dT_dz_prime=subs(dT_dz,[x,z],[x_prime,z_prime]);
%Define Green's function
Gxh= (1/(4*pi))*(3*(xm/(xm^2 + zp^2)) + 2*(xm*zm^2/(xm^2 +zm^2)^2))-(1/pi)*(3*(xm*(z_prime*zp + xm^2)/(xm^2 + zp^2)^2)-(3*(z_prime)^2*xm*zp*2 +xm^3*(4*z_prime^2 + 6*z*z_prime + z^2 + xm^2))/(xm^2+zp^2)^3);
Gxv= (-1/(4*pi))*((zp/(xm^2 + zp^2))+ 2*((xm^2*zm/(xm^2+zm^2)^2)-(xm^2*zm)/(xm^2 +zm^2)^2))-(1/(2*pi)*(2*(zp/(xm^2 + zp^2)))-((2*z-z_prime)*(zp^2-xm^2)/(xm^2+zp^2)^2)+(2*z*z_prime*zp*(3*xm^2-zp^2))/(xm^2+zp^2)^3);
p= P_l*exp(-(-v*t)/2*alpha)/(4*pi*k*(-v*t));
Gzh= (-1/(4*pi))*(3*(xm/(xm^2 +zp^2)- (xm/(xm^2 +zm^2))))- (1/pi)*(3*((xm*(z_prime*zp +xm^2)/(xm^2 +zp^2)^2))-(3*z_prime^2*xm*zp^2 + xm^3*(4*z_prime^2 + 6*z*z_prime +z^2 +xm^2)/(xm^2 +zp^2)^3));
Gxzh= (1/(4*pi))*((zp/(xm^2 + zp^2))-(zm/(xm^2 + zm^2))+2*((zp*xm^2/(xm^2 +zp^2)^2)-(zm*xm^2/(xm^2 +zm^2)^2)))-(1/pi)*(3*(z_prime*zp^2 +xm^2*(2*z+z_prime)/2*(xm^2+zp^2)^2)-(z_prime^3*(z_prime^2 + 3*z*z_prime + 3*z^2) + z^3*z_prime^2 + xm^2*(z_prime^3 + 6*z*z_prime^2 + 6*z^2*z_prime + z^3)+z*xm^4)/(xm^2 +zp^2)^3);
Gzv= (1/(4*pi))*(3*((zp/(xm^2 + zp^2))-(zm/(xm^2 + zm^2))) +2*((zm*xm^2/(xm^2 +zm^2)^2)-(zp*xm^2/(xm^2 + zp^2)^2))) - (1/(2*pi))*(2*zp/(xm^2+zp^2) + ((2*z+z_prime)*(zp^2-xm^2)/(xm^2 +zp^2)^2)-(2*z*z_prime*zp*(3*xm^2 -zp^2)/(xm^2 +zp^2)^3));
Gxzv= (xm/(4*pi))*((1/(xm^2+zp^2) - 1/(xm^2+zm^2)) + 2*((zp^2/(xm^2 +zp^2))-(zm^2/(xm^2+zm^2)^2)))-(xm/(2*pi))*((4*z*zp/(xm^2 +zp^2)^2) + (2*z*z_prime*zp*(3*zp^2 -xm^2)/(xm^2 + zp^2)^3));
%Convert to a function handle
T = matlabFunction(T);
p = matlabFunction(p);
fun1 = matlabFunction(Gxh * dT_dx + Gxv * dT_dz);
fun2 = matlabFunction(Gzh * dT_dx + Gzv * dT_dz);
fun3 = matlabFunction(Gxzh * dT_dx + Gxzv * dT_dz);
%Define terms as funciton handles for Sigma_xx_therm
term1_xx = @(t,x,z) integral2(@(x_prime,z_prime) fun1(t,x,z,x_prime,z_prime),-3,10,-3,0, 'AbsTol', 1e-6, 'RelTol', 1e-6);
term2_xx = @(x, z) (2 * z) / pi * integral(@(t) (p(t) .* (t - x).^2 ./ ((t - x).^2 + z^2).^2), 0, 0.1, 'AbsTol', 1e-6, 'RelTol', 1e-6);
term3_xx = @(t,x,z) -(alpha * E * T(t,x,z)) / (1 - 2*nu);
Sigma_xx_therm = -(alpha*E/(1-2*nu))*term1_xx(0,0,-0.0005) + term2_xx(0,-0.0001) + term3_xx(0,0,-0.0005);
%Define terms as funciton handles for Sigma_zz_therm
term1_zz = @(t,x,z) integral2(@(x_prime,z_prime) fun2(t,x,z,x_prime,z_prime),-3,10,-3,0, 'AbsTol', 1e-6, 'RelTol', 1e-6);
term2_zz = @(x, z) (2 * z^3) / pi * integral(@(t) (p(t) ./ ((t - x).^2 + z^2).^2), 0, 0.1, 'AbsTol', 1e-6, 'RelTol', 1e-6);
term3_zz = @(t,x,z) -(alpha * E * T(t,x,z)) / (1 - 2*nu);
Sigma_zz_therm = -(alpha*E/(1-2*nu))*term1_zz(0,0,-0.0005) + term2_zz(0,-0.0001) + term3_zz(0,0,-0.0005);
%Define terms as funciton handles for Sigma_xz_therm
term1_xz = @(t,x,z) integral2(@(x_prime,z_prime) fun3(t,x,z,x_prime,z_prime),-3,10,-3,0, 'AbsTol', 1e-6, 'RelTol', 1e-6);
term2_xz = @(x, z) (2 * z^2) / pi * integral(@(t) (p(t) .* (t - x)./ ((t - x).^2 + z^2).^2), 0, 0.1, 'AbsTol', 1e-6, 'RelTol', 1e-6);
term3_xz = @(t,x,z) -(alpha * E * T(t,x,z)) / (1 - 2*nu);
Sigma_xz_therm = -(alpha*E/(1-2*nu))*term1_xz(0,0,-0.0005) + term2_xz(0,-0.0001) + term3_xz(0,0,-0.0005);
% components of deviatoric stress
S_xx=(2*sigma_xx - sigma_yy - sigma_zz)/3;
S_yy= (2*sigma_yy - sigma_xx - sigma_zz)/3;
S_zz=(2*sigma_zz - sigma_yy - sigma_xx)/3;
S_eq=sqrt(sigma_xx^2 + sigma_yy^2 + sigma_zz^2);
n_xx=S_xx/S_eq;
n_yy=S_yy/S_eq;
n_zz=S_zz/S_eq;
%Equations
eqn1= (1/E)*(sigma_xx-nu*(sigma_yy+sigma_zz))+ (1/h)*(sigma_xx*n_xx + sigma_yy*n_yy + sigma_zz*n_zz)*n_xx == psi*((1/E)*(Sigma_xx_therm - nu*(sigma_yy + Sigma_zz_therm))+(1/h)*(Sigma_xx_therm*n_xx + sigma_yy*n_yy + Sigma_zz_therm*n_zz)*n_xx);
eqn2= (1/E)*(sigma_yy-nu*(sigma_xx+Sigma_zz_therm))+(1/h)*(sigma_xx*n_xx + sigma_yy*n_yy +sigma_zz*n_zz)*n_yy==0;
eqn3= sigma_yy==0.5*(sigma_xx+sigma_zz);
sol=solve([eqn1,eqn2,eqn3],[sigma_xx,sigma_yy,sigma_zz]);
sol1=sol.sigma_xx;
sol2=sol.sigma_yy;
sol3=sol.sigma_zz;
  댓글 수: 1
Dyuman Joshi
Dyuman Joshi 2023년 8월 28일
If you want to use solve, it would be better to stick to using symbolic expressions/functions, and avoid converting to function handles.

댓글을 달려면 로그인하십시오.

답변 (1개)

John D'Errico
John D'Errico 2023년 8월 28일
P_l=50;
v=0.1;
k=15;
Tm=1375;
T0=300;
alpha=3.75*10^(-6);
E= 190*10^9;
h= 100*10^9;
G= 150*10^9;
nu=0.3;
psi= 1- exp(-(0.45*h)/(2*G));
syms x z x_prime z_prime t dT_dx dT_dz sigma_xx sigma_yy sigma_zz;
xm=x-x_prime;
zp=z+z_prime;
zm=z-z_prime;
% Define the terms
T=(P_l*exp(-v*((sqrt((x-v*t)^2 + (z)^2) + (x-v*t)))/(2*alpha)))/(4*3.14*k*sqrt((x-v*t)^2 + (z)^2)) +T0;
dT_dx=diff(T,x);
dT_dx_prime=subs(dT_dx,[x,z],[x_prime,z_prime]);
dT_dz=diff(T,z);
dT_dz_prime=subs(dT_dz,[x,z],[x_prime,z_prime]);
%Define Green's function
Gxh= (1/(4*pi))*(3*(xm/(xm^2 + zp^2)) + 2*(xm*zm^2/(xm^2 +zm^2)^2))-(1/pi)*(3*(xm*(z_prime*zp + xm^2)/(xm^2 + zp^2)^2)-(3*(z_prime)^2*xm*zp*2 +xm^3*(4*z_prime^2 + 6*z*z_prime + z^2 + xm^2))/(xm^2+zp^2)^3);
Gxv= (-1/(4*pi))*((zp/(xm^2 + zp^2))+ 2*((xm^2*zm/(xm^2+zm^2)^2)-(xm^2*zm)/(xm^2 +zm^2)^2))-(1/(2*pi)*(2*(zp/(xm^2 + zp^2)))-((2*z-z_prime)*(zp^2-xm^2)/(xm^2+zp^2)^2)+(2*z*z_prime*zp*(3*xm^2-zp^2))/(xm^2+zp^2)^3);
p= P_l*exp(-(-v*t)/2*alpha)/(4*pi*k*(-v*t));
Gzh= (-1/(4*pi))*(3*(xm/(xm^2 +zp^2)- (xm/(xm^2 +zm^2))))- (1/pi)*(3*((xm*(z_prime*zp +xm^2)/(xm^2 +zp^2)^2))-(3*z_prime^2*xm*zp^2 + xm^3*(4*z_prime^2 + 6*z*z_prime +z^2 +xm^2)/(xm^2 +zp^2)^3));
Gxzh= (1/(4*pi))*((zp/(xm^2 + zp^2))-(zm/(xm^2 + zm^2))+2*((zp*xm^2/(xm^2 +zp^2)^2)-(zm*xm^2/(xm^2 +zm^2)^2)))-(1/pi)*(3*(z_prime*zp^2 +xm^2*(2*z+z_prime)/2*(xm^2+zp^2)^2)-(z_prime^3*(z_prime^2 + 3*z*z_prime + 3*z^2) + z^3*z_prime^2 + xm^2*(z_prime^3 + 6*z*z_prime^2 + 6*z^2*z_prime + z^3)+z*xm^4)/(xm^2 +zp^2)^3);
Gzv= (1/(4*pi))*(3*((zp/(xm^2 + zp^2))-(zm/(xm^2 + zm^2))) +2*((zm*xm^2/(xm^2 +zm^2)^2)-(zp*xm^2/(xm^2 + zp^2)^2))) - (1/(2*pi))*(2*zp/(xm^2+zp^2) + ((2*z+z_prime)*(zp^2-xm^2)/(xm^2 +zp^2)^2)-(2*z*z_prime*zp*(3*xm^2 -zp^2)/(xm^2 +zp^2)^3));
Gxzv= (xm/(4*pi))*((1/(xm^2+zp^2) - 1/(xm^2+zm^2)) + 2*((zp^2/(xm^2 +zp^2))-(zm^2/(xm^2+zm^2)^2)))-(xm/(2*pi))*((4*z*zp/(xm^2 +zp^2)^2) + (2*z*z_prime*zp*(3*zp^2 -xm^2)/(xm^2 + zp^2)^3));
%Convert to a function handle
T = matlabFunction(T);
p = matlabFunction(p);
fun1 = matlabFunction(Gxh * dT_dx + Gxv * dT_dz);
fun2 = matlabFunction(Gzh * dT_dx + Gzv * dT_dz);
fun3 = matlabFunction(Gxzh * dT_dx + Gxzv * dT_dz);
%Define terms as funciton handles for Sigma_xx_therm
term1_xx = @(t,x,z) integral2(@(x_prime,z_prime) fun1(t,x,z,x_prime,z_prime),-3,10,-3,0, 'AbsTol', 1e-6, 'RelTol', 1e-6);
term2_xx = @(x, z) (2 * z) / pi * integral(@(t) (p(t) .* (t - x).^2 ./ ((t - x).^2 + z^2).^2), 0, 0.1, 'AbsTol', 1e-6, 'RelTol', 1e-6);
term3_xx = @(t,x,z) -(alpha * E * T(t,x,z)) / (1 - 2*nu);
Sigma_xx_therm = -(alpha*E/(1-2*nu))*term1_xx(0,0,-0.0005) + term2_xx(0,-0.0001) + term3_xx(0,0,-0.0005);
%Define terms as funciton handles for Sigma_zz_therm
term1_zz = @(t,x,z) integral2(@(x_prime,z_prime) fun2(t,x,z,x_prime,z_prime),-3,10,-3,0, 'AbsTol', 1e-6, 'RelTol', 1e-6);
term2_zz = @(x, z) (2 * z^3) / pi * integral(@(t) (p(t) ./ ((t - x).^2 + z^2).^2), 0, 0.1, 'AbsTol', 1e-6, 'RelTol', 1e-6);
term3_zz = @(t,x,z) -(alpha * E * T(t,x,z)) / (1 - 2*nu);
Sigma_zz_therm = -(alpha*E/(1-2*nu))*term1_zz(0,0,-0.0005) + term2_zz(0,-0.0001) + term3_zz(0,0,-0.0005);
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 9.8e+16. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.
%Define terms as funciton handles for Sigma_xz_therm
term1_xz = @(t,x,z) integral2(@(x_prime,z_prime) fun3(t,x,z,x_prime,z_prime),-3,10,-3,0, 'AbsTol', 1e-6, 'RelTol', 1e-6);
term2_xz = @(x, z) (2 * z^2) / pi * integral(@(t) (p(t) .* (t - x)./ ((t - x).^2 + z^2).^2), 0, 0.1, 'AbsTol', 1e-6, 'RelTol', 1e-6);
term3_xz = @(t,x,z) -(alpha * E * T(t,x,z)) / (1 - 2*nu);
Sigma_xz_therm = -(alpha*E/(1-2*nu))*term1_xz(0,0,-0.0005) + term2_xz(0,-0.0001) + term3_xz(0,0,-0.0005);
% components of deviatoric stress
S_xx=(2*sigma_xx - sigma_yy - sigma_zz)/3;
S_yy= (2*sigma_yy - sigma_xx - sigma_zz)/3;
S_zz=(2*sigma_zz - sigma_yy - sigma_xx)/3;
S_eq=sqrt(sigma_xx^2 + sigma_yy^2 + sigma_zz^2);
n_xx=S_xx/S_eq;
n_yy=S_yy/S_eq;
n_zz=S_zz/S_eq;
%Equations
eqn1= (1/E)*(sigma_xx-nu*(sigma_yy+sigma_zz))+ (1/h)*(sigma_xx*n_xx + sigma_yy*n_yy + sigma_zz*n_zz)*n_xx == psi*((1/E)*(Sigma_xx_therm - nu*(sigma_yy + Sigma_zz_therm))+(1/h)*(Sigma_xx_therm*n_xx + sigma_yy*n_yy + Sigma_zz_therm*n_zz)*n_xx);
eqn2= (1/E)*(sigma_yy-nu*(sigma_xx+Sigma_zz_therm))+(1/h)*(sigma_xx*n_xx + sigma_yy*n_yy +sigma_zz*n_zz)*n_yy==0;
eqn3= sigma_yy==0.5*(sigma_xx+sigma_zz);
sol=solve([eqn1,eqn2,eqn3],[sigma_xx,sigma_yy,sigma_zz]);
sol.sigma_xx
ans = 
vpa(sol.sigma_xx)
ans = 
So it looks like three roots, Two of which are complex conjugates. Your problem is probably implicitly a cubic polynomial, so that should not be a surprise.
vpa(sol.sigma_yy)
ans = 
vpa(sol.sigma_zz)
ans = 
If you want only the real roots, then keep only root number 1. Where is the problem?
As far as plotting a solution as a function of z, since z is not present in the solution, that does not make complete sense.

카테고리

Help CenterFile Exchange에서 Programming에 대해 자세히 알아보기

제품


릴리스

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by