This code is correct theoratically, but cannot be runned.

1 次查看（过去 30 天）
Zhe Liu
Zhe Liu 2018-8-19
评论： Stephen23 2018-8-20
Hello, I tried to run this code to get a 3d graph of a deformation of a fiber, but I get a lot of NaN, when I run it. I checked but did not find /0. Could someone help me to check and fix it? Thanks!!!
Just to be clear...This code is the result of solving a pde and get the deformation shape of a fiber.
Here is my code:
if true
%
clear
clc
m0 = 5;
l = pi;
c = 1;
alpha = 150;
[x,y] = meshgrid(-.5:.1:.5,-.5:.1:.5);
dimx = size(x);
dimy = size(y);
lengthx = dimx(1,1);
lengthy = dimy(1,2);
u1 = zeros(lengthx,lengthy);
u2 = zeros(lengthx,lengthy);
for m = 1:2:5
pa1 = sqrt(-((sqrt(1 - 4.*alpha.*m.^2) - 1)./(2.*alpha)));
pa2 = sqrt((sqrt(1 - 4.*alpha.*m.^2) + 1)./(2.*alpha));
%kvalue
k11 = sin(m.*y).*(pa2.^3.*exp(pa2.*x) + pa2.^3.*exp(-pa2.*x));
k12 = sin(m.*y).*(pa1.^3.*exp(pa1.*x) + pa1.^3.*exp(-pa1.*x));
k21 = -m.*cos(m.*y).*(pa2.^2.*exp(c.*pa2) - pa2.^2.*exp(-c.*pa2)) ;
k22 = -m.*cos(m.*y).*(pa1.^2.*exp(c.*pa1) - pa1.^2.*exp(-c.*pa1)) ;
detK = (k11*k22 - k12*k21);
%fs
A = (4*m0/(m*pi));
B = (-1)^(m/2 - 0.5);
C = m;
fs = A*B*cos((C * x));
%unknowns
C1 = -(fs*k12)/detK;
C3 = (fs*k11)/detK;
C2 = - C1;
C4 = - C3;
%u1u2
u1 = u1 + m.*sin(m.*y).*(C1.*exp(pa2.*x) + C3.*exp(pa1.*x) + C2.*exp(-pa2.*x) + C4.*exp(-pa1.*x));
u2 = u2 + cos(m.*y).*(C1.*pa2.*exp(pa2.*x) + C3.*pa1.*exp(pa1.*x) - C2.*pa2.*exp(-pa2.*x) - C4.*pa1.*exp(-pa1.*x));
end
figure
plot(y,u1(:,1))
title('In-plane deflection of the membrane: u1(x,y)')
xlabel('x coordinate (nm)')
ylabel('y along the cross section x=0')
%grid on
figure
plot(y,u2(:,1))
title('In-plane deflection of the membrane: u2(x,y)')
xlabel('x coordinate (nm)')
ylabel('y along the cross section x=0')
%grid on
% plot contour surface
figure;
surf(x,y,u1)
shading('interp')
colorbar
%%colormap gray
title('Membrane in-plane Deflection');
xlabel('x coordinate (nm)');
ylabel('y coordinate (nm)');
zlabel('u1(x,y) (nm)');
% plot contour surface
figure;
surf(x,y,u2)
shading('interp')
colorbar
%%colormap gray
title('Membrane in-plane Deflection');
xlabel('x coordinate (nm)');
ylabel('y coordinate (nm)');
zlabel('u2(x,y) (nm)');
end

请先登录，再回答此问题。

回答（1 个）

Image Analyst
Image Analyst 2018-8-19
Your pa1 and pa2 are complex numbers. And your det is an array, not a single numbers. I suggest you take a closer look at the actual numbers and formulas.
  1 个评论
Stephen23
Stephen23 2018-8-20
Zhe Liu's "Flag" moved here:
Thank you very much! Pa1 and pa2 come from the solution of a pde. I would double check it.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Visualization 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by