上机实验2:五点差分格式法 偏微分方程(Matlab)实验报告
——五点差分格式法 一、 实验题目
设G是形如下图的十字形域,由五个相等的单位正方形组成,用五点差分格式求下列边值问题的数值解:
?u?u?2??1,于G2?x?yu=0,于?G二、 实验原理
22G1 G 取定沿X轴和Y轴方向的步长h1和h2,h??h?h221122?,作两族与坐
标轴平行的直线:x=ih1,y=jh2,i,j?0,?1,?2,?
若(xi,yj)为正则内点,沿x,y方向分别用二阶中心差商代替
uxx和uyy则得
?[ui?1,j?2uij?ui?1,jh12?ui,j?1?2uij?ui,j?12h2]?fij
特别取正方形网格:h1?h2?h,则原差分方程可简化为
1h2uij?(ui?1,j?ui,j?1?ui?1,j?ui,j?1)?fij
44三、 实验程序
1)function uxy = EllIni2Uxl(x,y)
format long;
uxy = 0;
2)function uxy = EllIni2Uxr(x,y) format long; uxy = y*(2-y);
3)function uxy = EllIni2Uyl(x,y) format long; uxy = 0;
4)function uxy = EllIni2Uyr(x,y) format long; if x < 1
uxy = x; else
uxy = 2 - x; end
5)function u = peEllip5(nx,minx,maxx,ny,miny,maxy) format long;
hx = (maxx-minx)/(nx-1); hy = (maxy-miny)/(ny-1); u0 = zeros(nx,ny); for j=1:ny
u0(j,1) = EllIni2Uxl(minx,miny+(j-1)*hy); u0(j,nx) = EllIni2Uxr(maxx,miny+(j-1)*hy); end
for j=1:nx
u0(1,j) = EllIni2Uyl(minx+(j-1)*hx,miny); u0(ny,j) = EllIni2Uyr(minx+(j-1)*hx,maxy); end
A = -4*eye((nx-2)*(ny-2),(nx-2)*(ny-2)); b = ones((nx-2)*(ny-2),1).*(-1); for i=1:(nx-2)*(ny-2) if mod(i,nx-2) == 1 if i==1
A(1,2) = 1; A(1,nx-1) = 1;
b(1) = - u0(1,2) - u0(2,1); else
if i == (ny-3)*(nx-2)+1 A(i,i+1) = 1; A(i,i-nx+2) = 1;
b(i) = - u0(ny-1,1) - u0(ny,2); else
A(i,i+1) = 1; A(i,i-nx+2) = 1; A(i,i+nx-2) = 1;
b(i) = - u0(floor(i/(nx-2))+2,1); end end else
if mod(i,nx-2) == 0 if i == nx-2
A(i,i-1) = 1; A(i,i+nx-2) = 1;
b(i) = - u0(1,nx-1) - u0(2,nx); else
if i == (ny-2)*(nx-2) A(i,i-1) = 1; A(i,i-nx+2) = 1;
b(i) = - u0(ny-1,nx) - u0(ny,nx-1); else
A(i,i-1) = 1; A(i,i-nx+2) = 1; A(i,i+nx-2) = 1;
b(i) = - u0(floor(i/(nx-2))+1,nx); end end else
if i>1 && i< nx-2 A(i,i-1) = 1; A(i,i+nx-2) = 1; A(i,i+1) = 1;
b(i) = - u0(1,i+1); else
if i > (ny-3)*(nx-2) && i < (ny-2)*(nx-2) A(i,i-1) = 1; A(i,i-nx+2) = 1; A(i,i+1) = 1;
b(i) = - u0(ny,mod(i,(nx-2))+1); else
A(i,i-1) = 1; A(i,i+1) = 1; A(i,i+nx-2) = 1; A(i,i-nx+2) = 1; end
end end end end
ul = A\\b;
for i=1:(ny-2) for j=1:(nx-2)
u(i,j) = ul((i-1)*(nx-2)+j); end end
format short;
四、 实验结果
>> u=peEllip5(25,0,3,25,0,3) u =
Columns 1 through 6
1.1448 2.2896 3.2671 2.2896 4.7466 6.7036 3.2670 6.7035 9.5050 4.0749 8.2955 11.8064 4.7370 9.5973 13.6990 5.2759 10.6576 15.2465 5.7089 11.5108 16.4951 6.0490 12.1816 17.4787 6.3055 12.6880 18.2221 6.4851 13.0426 18.7432 6.5922 13.2542 19.0543 6.6294 13.3278 19.1623 6.5977 13.2651 19.0700 6.4964 13.0648 18.7754 6.3231 12.7225 18.2720 6.0736 12.2298 17.5486 5.7415 11.5747 16.5879 5.3175 10.7395 15.3659 4.7891 9.7000 13.8492 4.1389 8.4221 11.9923 3.3445 6.8572 9.7318 2.3820 4.9303 6.9762 1.2532 2.5058 3.5896
Columns 7 through 12
4.0750 8.2957 11.8066 14.7262 17.1457 19.1343 20.7444 22.0159 22.9785 23.6540 24.0573 24.1973 24.0771 23.6946 23.0417 22.1047 20.8626 19.2868 17.3385 14.9662 12.1014 8.6530 4.5016 4.7372 9.5977 13.6994 17.1460 20.0233 22.4005 24.3324 25.8619 27.0220 27.8369 28.3237 28.4924 28.3467 27.8841 27.0957 25.9657 24.4710 22.5803 20.2519 17.4326 14.0546 10.0330 5.2637 5.2763 10.6583 15.2474 19.1351 22.4011 25.1120 27.3228 29.0775 30.4106 31.3481 31.9082 32.1019 31.9332 31.3995 30.4911 29.1913 27.4756 25.3113 22.6562 19.4577 15.6515 11.1607 5.8950
5.7095 6.0499 6.3068 6.4869 6.5946 6.6327 11.5119 12.1833 12.6905 13.0461 13.2590 13.3341 16.4966 17.4810 18.2256 18.7482 19.0610 19.1713 20.7461 22.0186 22.9826 23.6600 24.0656 24.2085 24.3339 25.8647 27.0264 27.8435 28.3330 28.5051 27.3238 29.0798 30.4148 31.3547 31.9178 32.1153 29.7693 31.7160 33.1982 34.2426 34.8684 35.0871 31.7146 33.8168 35.4193 36.5492 37.2260 37.4615 33.1947 35.4172 37.1129 34.2366 36.5445 38.3062 34.8591 37.2178 39.0186 35.0737 37.4491 39.2621 34.8848 37.2427 39.0411 34.2896 36.5959 38.3529 33.2780 35.4985 37.1872 31.8328 33.9329 35.5265 29.9289 31.8739 33.3454 27.5331 29.2885 30.6112 24.6039 26.1358 27.2838 21.0904 22.3671 23.3156 16.9329 17.9266 18.6526 12.0633 12.7537 13.2351 6.4058 6.7898 6.9996
Columns 13 through 18
6.6020 6.5019 6.3302 13.2733 13.0755 12.7361 19.0818 18.7906 18.2915 24.0918 23.7135 23.0660 28.3635 27.9059 27.1234 31.9512 31.4230 30.5211 34.9031 34.3137 33.3090 37.2603 36.6196 35.5294 39.0571 38.3751 37.2167 40.3201 39.6070 38.3981 41.0688 40.3347 39.0928 41.3152 40.5704 39.3129 41.0640 40.3185 39.0632 40.3124 39.5766 38.3409 39.0503 38.3344 37.1363 37.2600 36.5745 35.4318 34.9165 34.2716 33.2024
38.3090 39.5489 40.2884 40.5395 40.3067 39.5873 38.3711 36.6403 34.3698 31.5272 28.0725 23.9590 19.1329 13.5347 7.0985 6.0826 12.2473 17.5734 22.1354 26.0007 29.2291 31.8719 33.9721 35.5644 36.6757 37.3254 37.5254 37.2803 36.5877 35.4380 33.8140 31.6909 39.0249 40.2920 41.0465 41.3007 41.0590 40.3187 39.0695 37.2938 34.9665 32.0552 28.5202 24.3148 19.3854 13.6723 7.1096 5.7530 11.5970 16.6195 20.9014 24.5150 27.5227 29.9775 31.9226 33.3930 34.4151 35.0078 35.1829 34.9450 34.2917 33.2139 31.6951 29.7112 39.2728 40.5475 41.3051 41.5577 41.3099 40.5588 39.2944 37.4988 35.1472 32.2070 28.6384 24.3945 19.4216 13.6596 7.0427 5.3325 10.7682 15.4061 19.3358 22.6350 25.3694 27.5926 29.3480 30.6698 31.5837 32.1078 32.2534 32.0249 31.4203 30.4309 29.0414 27.2288
31.9871 31.3931 30.4153 29.0362 27.2300 24.9623 28.4318 27.8985 27.0295 25.8084 24.2105 22.2018 24.2032 23.7396 22.9959 21.9575 20.6016 18.8966 19.2468 18.8607 18.2571 17.4239 16.3420 14.9838 13.5018 13.1993 12.7480 12.1390 11.3586 10.3859 6.9016 6.6868 6.3964 6.0257 5.5674 5.0103
Columns 19 through 23
4.8087 9.7373 13.9009 17.4005 20.3200 22.7271 24.6754 26.2071 27.3546 28.1421 28.5864 28.6980 28.4809 27.9336 27.0481 25.8106 24.2003 22.1885 19.7377 16.7993 13.3107 9.1911 4.3378
>>
4.1651 8.4714 12.0597 15.0454 17.5175 19.5436 21.1749 22.4504 23.3995 24.0436 24.3978 24.4711 24.2673 23.7849 23.0174 21.9528 20.5732 18.8537 16.7613 14.2523 11.2685 7.7299 3.5247 3.3801 6.9236 9.8209 12.2038 18.7492 18.9214 13.2014 8.7810 5.9352 2.5313 2.4318 5.0218 7.0967 8.7880 10.1681 11.2845 12.1707 12.8507 13.3421 13.6577 13.8060 13.7923 13.6187 13.2846 12.7864 12.1176 11.2686 10.2260 8.9724 7.4861 5.7401 3.6985 1.2905 1.3254 2.6352 3.7561 4.6832 5.4388 6.0445 6.5173 6.8697 7.1109 7.2475 7.2840 7.2230 7.0657 6.8118 6.4595 6.0055 5.4449 4.7711 3.9762 3.0518 1.9949 0.8282 -0.3180 14.1611 15.7548 17.0301 18.0201 19.2352 19.4901 19.5214 19.3322 18.2838 17.4100 16.2859 14.8917 11.1802