Am I using poicalc correctly?

조회 수: 1 (최근 30일)
William Kett
William Kett 2021년 10월 4일
편집: William Kett 2021년 10월 4일
The solution returned using poicalc is too large by almost 1 order or magnitude. I expect the maximum z value of the surface to be 1, but instead I am getting 7.9976 (almost 8). I will post exactly what I did below.
I used the periodic, 2D function
And I set a = 1 to make things simpler
syms u(x,y,a)
u(x,y,a) = 2^(4*a) * x^a * (1-x)^a * y^a * (1-y)^a;
u(x,y,a) = subs(u(x,y,a),a,1); %sub in a = 1 into equation
u(x,y,a)
ans = 
%% Cartesian Space
% x = [0,.1,.2,...,X=1], y = [0,.1,.2,...,Y]
X = 1; Y = 1;
Lx = 99; Ly = 99;
dx = X/Lx; dy = Y/Ly;
x_arr = linspace(0,X,Lx+1);
y_arr = transpose(linspace(0,Y,Ly+1));
x_arr
x_arr = 1×100
0 0.0101 0.0202 0.0303 0.0404 0.0505 0.0606 0.0707 0.0808 0.0909 0.1010 0.1111 0.1212 0.1313 0.1414 0.1515 0.1616 0.1717 0.1818 0.1919 0.2020 0.2121 0.2222 0.2323 0.2424 0.2525 0.2626 0.2727 0.2828 0.2929
y_arr
y_arr = 100×1
0 0.0101 0.0202 0.0303 0.0404 0.0505 0.0606 0.0707 0.0808 0.0909
u_func(x,y,a) = u(x,y,a); % This is a term that we will be using to sub into u(x,y,a)
u(x,y,a) = subs(u(x,y,a),x,x_arr(1,:)); % sub in x
u(x,y,a) = subs(u(x,y,a),y,y_arr(:,1)); % sub in y. I noticed that I could sub in a transposed array instead of using a meshgrid
u_xy = double(u(x,y,a)) % We need to convert this to a double for plotting purposes later on
u_xy = 100×100
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0016 0.0032 0.0047 0.0062 0.0077 0.0091 0.0105 0.0119 0.0132 0.0145 0.0158 0.0170 0.0182 0.0194 0.0206 0.0217 0.0228 0.0238 0.0248 0.0258 0.0267 0.0277 0.0285 0.0294 0.0302 0.0310 0.0317 0.0325 0.0331 0 0.0032 0.0063 0.0093 0.0123 0.0152 0.0180 0.0208 0.0235 0.0262 0.0288 0.0313 0.0337 0.0361 0.0385 0.0407 0.0429 0.0450 0.0471 0.0491 0.0511 0.0529 0.0547 0.0565 0.0582 0.0598 0.0613 0.0628 0.0642 0.0656 0 0.0047 0.0093 0.0138 0.0182 0.0225 0.0268 0.0309 0.0349 0.0389 0.0427 0.0464 0.0501 0.0536 0.0571 0.0604 0.0637 0.0669 0.0699 0.0729 0.0758 0.0786 0.0813 0.0839 0.0863 0.0887 0.0910 0.0933 0.0954 0.0974 0 0.0062 0.0123 0.0182 0.0241 0.0297 0.0353 0.0408 0.0461 0.0513 0.0563 0.0613 0.0661 0.0708 0.0753 0.0798 0.0841 0.0882 0.0923 0.0962 0.1000 0.1037 0.1072 0.1106 0.1139 0.1171 0.1201 0.1230 0.1258 0.1285 0 0.0077 0.0152 0.0225 0.0297 0.0368 0.0437 0.0504 0.0570 0.0634 0.0697 0.0758 0.0817 0.0875 0.0932 0.0986 0.1040 0.1091 0.1141 0.1190 0.1237 0.1282 0.1326 0.1368 0.1409 0.1448 0.1486 0.1522 0.1556 0.1589 0 0.0091 0.0180 0.0268 0.0353 0.0437 0.0519 0.0599 0.0677 0.0753 0.0827 0.0900 0.0970 0.1039 0.1106 0.1171 0.1234 0.1296 0.1355 0.1413 0.1468 0.1522 0.1574 0.1625 0.1673 0.1719 0.1764 0.1807 0.1848 0.1887 0 0.0105 0.0208 0.0309 0.0408 0.0504 0.0599 0.0691 0.0781 0.0869 0.0955 0.1038 0.1120 0.1199 0.1276 0.1352 0.1425 0.1495 0.1564 0.1630 0.1695 0.1757 0.1817 0.1875 0.1931 0.1984 0.2036 0.2085 0.2132 0.2178 0 0.0119 0.0235 0.0349 0.0461 0.0570 0.0677 0.0781 0.0883 0.0982 0.1079 0.1174 0.1266 0.1356 0.1443 0.1528 0.1610 0.1690 0.1768 0.1843 0.1916 0.1986 0.2054 0.2120 0.2183 0.2243 0.2301 0.2357 0.2411 0.2462 0 0.0132 0.0262 0.0389 0.0513 0.0634 0.0753 0.0869 0.0982 0.1093 0.1201 0.1306 0.1409 0.1508 0.1606 0.1700 0.1792 0.1881 0.1967 0.2051 0.2132 0.2210 0.2285 0.2358 0.2428 0.2496 0.2561 0.2623 0.2682 0.2739
%% Now we plot the solution
figure(1)
contour3(x_arr,y_arr,u_xy);
surface(x_arr,y_arr,u_xy);
title('u(x,y)');
xlabel('x = linspace(0,1,11)');ylabel('y = linspace(0,1,11)');
Then we set up the Poisson equation that we are trying to solve for in order to get back to the solution graphed above.
f_func = diff(u_func,x,2) + diff(u_func,y,2);
A = subs(f_func,'x',x_arr);
B = subs(A,'y', y_arr);
f = double(B);
size(f)
ans = 1×2
100 100
Now that we have f(x,y)< we will use poicalc on f(x,y) to get back to u(x,y).
%% Use poicalc function
% Now we run the poicalc function which takes in parameters (f,h1,h2,n1,n2)
h1 = dx; % gridspace_x
h2 = dy; % gridspace_y
n1 = 100; %number of points in x. supposed to be 100 but n1*n2 has to equal the number of rows in f
n2 = 1; % number of points in y. Supposed to be 100.
% n2 seems to dictate the number of local maxima and minima or
% slopes
% if I put n2 = 2, I get 2 slopes, and if n1 = n2 = 10, I get 10
% slopes.
u_calc = poicalc(f,h1,h2,n1,n2);
%% Now we plot it to examine how well it matches our first plot
figure(2)
contour3(y_arr,x_arr,abs(u_calc));
surface(y_arr,x_arr,abs(u_calc));
title('u(x,y) computed using poicalc function');
xlabel('x');ylabel('y');
max(max(abs(u_xy)))
ans = 0.9998
max(max(abs(u_calc)))
ans = 7.9976
Personally I don't have too much issue with the shape, my biggest issue is with the magnitude being off. The maximum value for u_xy in the first plot is very nearly 1, and the max for the solution obtained for u(x,y) here is very nearly 8. What exactly am I doing wrong with the poicalc function that is resulting in this disparity in magnitude?

답변 (0개)

카테고리

Help CenterFile Exchange에서 Operators and Elementary Operations에 대해 자세히 알아보기

제품


릴리스

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by