最近在学习FLAC软件，读到几篇文章，感觉这几篇文章在数值模拟学习中具有指导性意义。特附上文章，仅供大家参考，交流。 1. Towards a methodology for rock mechanics modeling.pdf （主要讲述在岩石力学数值中如何正确看待和使用数值模拟软件） 2. VerificationValidation and Confirmation of Numerical models in the earth sciences.pdf （主要探讨了在地球科学数值模拟中3个关键但滥用的英文单词；兼谈如何正确看待数值模拟软件） 3. 试验结果的数学拟合与力学模型.pdf （主要探讨了如何合理使用岩石力学的实验数据，结合实例解释了“过度拟合必将失真”） 4. Numercial model calibration for simulating coal pillars, gob and overburden response.pdf （结合实例，探讨了如何在采矿数值模拟中选取合理的数值，及如何对参数进行校正） 5. Geomech Models - Keep simple.pdf （采矿中简化数值模型的意见和价值）
I am an ocean modeler by training. I know what my model can and cannot do. I try my best to validate my model output using observations, so should you.
需要上机练习编程：差分法解拉普拉斯方程的第一边值问题。 自己编制的程序如下： 文件名：Lap-Eq Numerical answer.m clc;clear; tic N=50 %划分的网格数====================== for m=1:N n=1:N-1; u(m,n)=0; u(m,N)=sin((m-1)*pi/(N-1)); end %定义边界条件======================= delta=ones(N,N); while delta1e-6 for m=2:N-1 n=2:N-1; a(m,n)=u(m,n); u(m,n)=(u(m+1,n)+u(m-1,n)+u(m,n+1)+u(m,n-1))/4; delta(m,n)=abs(u(m,n)-a(m,n))/u(m,n); end end X=1:N;Y=1:N; mesh(X,Y,u(X,Y)) toc 所用的计算时间为Elapsed time is 3.672000 seconds. 1.考虑程序中的循环控制条件while delta=10e-6的意义。 经过单步调试，得知这个表达式只是对最后一个delta进行比较，而不是所有的delta，因此并不满足计算条件。结果是错误的。要求每个计算点的delta都要10e-6，因此需要该在程序。 clc;clear; tic N=50 %划分的网格数====================== for m=1:N n=1:N-1; u(m,n)=0; u(m,N)=sin((m-1)*pi/(N-1)); end %定义边界条件======================= delta=ones(N,N); for m=2:N-1 n=2:N-1; while delta(m,n)1e-4 for m=2:N-1 n=2:N-1; a(m,n)=u(m,n); u(m,n)=(u(m+1,n)+u(m-1,n)+u(m,n+1)+u(m,n-1))/4; delta(m,n)=abs(u(m,n)-a(m,n))/u(m,n); end end end X=1:N;Y=1:N; mesh(X,Y,u(X,Y)) delta delta1e-6 toc 这样一来，所有的点都满足了。但是这种算法做了太多的冗余计算。对每个点的delta分别调到误差范围，所作的计算次数太多太多了。从时间可看出 Elapsed time is 144.610000 seconds.需要进行算法改进，考虑每次计算结束得到一个delta的矩阵，只要矩阵中的最大者满足误差范围则所有的点都满足了，因此改为： clc;clear; tic N=50 %划分的网格数====================== for m=1:N n=1:N-1; u(m,n)=0; u(m,N)=sin((m-1)*pi/(N-1)); end %定义边界条件======================= delta=zeros(N,N); maxd=1; while maxd1e-4 for m=2:N-1 n=2:N-1; a(m,n)=u(m,n); u(m,n)=(u(m+1,n)+u(m-1,n)+u(m,n+1)+u(m,n-1))/4; delta(m,n)=abs(u(m,n)-a(m,n))/u(m,n); end maxd=max(delta(:)); end X=1:N;Y=1:N; mesh(X,Y,u(X,Y)) maxd1e-4 toc 现在就完全解决了上述问题了。Elapsed time is 1.954000 seconds. 若误差要求改为1e-5，则运行时间为Elapsed time is 3.797000 seconds. 若误差要求改为1e-6，则运行时间为Elapsed time is 5.687000 seconds. 若划分的网格节点数N=500， tolerance=1e-5= Elapsed time is 3794.234000 seconds. 几点需要说明的： 1) 对于二维或多维矩阵，找其最大值的表达式为max(A(:))，A代表矩阵名称。A(:)代表矩阵A的所有元素以单序号方式引用。这样找到的最大值才是一个数值。若单纯的使用max(A)则对于二维矩阵会得到一个行矩阵，对应于A中每列的最大值。 2) maxd1e-4的作用，是检验是否所有的delta都已经满足误差要求了。若满足，该式子的返回值为1，即为真。 另外，改变delta的存储情况也可减少存储空间，加快计算。以下是不用矩阵存储delta， 因为我们不需要知道每个delta值的表现形式，因此可以对每个delta进行比较只用一个值来存储它。程序如下： %椭圆型方程的数值计算典型例题，Laplace方程的第一边值问题。 clc;clear;clf; tic N=5 %划分的网格节点数 tol=1e-5 %差分误差tolerance要求 %计算精度控制参量====================== for m=1:Nn=1:N-1; u(m,n)=0; u(m,N)=sin((m-1)*pi/(N-1)); end %定义边界条件======================= delta=1;% 用于存储两次计算的相对误差 %maxd=1;%N*N 个相对误差中最大的一个 toi=0;%times of iteration迭代计算次数 while deltatol for m=2:N-1 n=2:N-1; a(m,n)=u(m,n); u(m,n)=(u(m+1,n)+u(m-1,n)+u(m,n+1)+u(m,n-1))/4; delta=max(abs(u(m,n)-a(m,n))/u(m,n)); % 注意这种形式的意义。 end toi=toi+1; %maxd=max(delta(:)); end X=1:N;Y=1:N; u mesh(X,Y,u(X,Y)) deltatol toi toc 另外，在进行判断时对delta进行的比较也可这样编写： while deltatol delta=0; for m=2:N-1 n=2:N-1; a(m,n)=u(m,n); u(m,n)=(u(m+1,n)+u(m-1,n)+u(m,n+1)+u(m,n-1))/4; delta=max(abs(u(m,n)-a(m,n))/u(m,n),delta); end toi=toi+1; %maxd=max(delta(:)); end 附：找矩阵中的最大值及其位置。 分情况（一维、二维或三维）而言： i.一维阵， =max(A)即可，a为最大值，b为位置; ii.二维矩阵 a=max(A(:)); =find(A==a); iii.三维 a=max(A(:)); ind = find(a==max(a(:))); = ind2sub( ,ind);
标签: numerical
