GPT辅助科学计算编程：两个实例拆解提示词设计与验证-北京尧图网络

从去年开始我做材料计算方向的Python脚本基本都在GPT辅助下完成这个系列也写到了第四篇。前面的内容讲了不少提示词框架和基础技巧今天这篇我打算完全换一种讲法直接拿两个计算力学和材料计算里的典型任务庖丁解牛一样把提示词怎么设计、代码怎么演进、验证怎么做、坑怎么踩全过程摊开来看。先说说我为什么坚持在计算材料这个方向上用GPT。早期我也用过它写通用程序说实话体验一般经常给出看起来能跑但物理上一塌糊涂的代码。直到我试了一次让它帮我推层合板刚度矩阵的公式并转成程序那一次虽然也出了错但错的位置非常有启发——它把Q矩阵里某个除法理解错了导致结果整体差一个倍数。从那时起我就意识到科学计算场景下GPT的正确用法不是我提需求它交代码而是我提供足够的物理约束和验证锚点让它在一个不太容易跑偏的窄轨道里工作。这个窄轨道就是提示工程要解决的问题。今天这篇我准备了两个实例一个2D桁架有限元求解器一个Lennard-Jones势函数下的原子能量与受力计算。一个偏力学一个偏材料正好覆盖标题里计算材料科学与力学两大块。每个实例我都会给出实际用过的提示词、跑出来的代码、以及中间踩过的坑希望能给正在用GPT辅助写科学计算代码的朋友一点参考。1. 科学计算的提示词为什么不能照搬通用提问法1.1 通用提问与科学计算提问的本质区别很多人用GPT写科学计算代码时习惯还是通用开发那套——用Python写一个计算梁弯曲变形的函数然后等着看结果。这种问法放在通用业务代码上勉强可以放在科学计算上几乎是必然翻车因为科学计算代码的正确性根本不取决于语法而取决于物理模型、数值格式、边界条件、单位制这些东西有没有被准确描述。我做过一个对比实验。第一个问题我让GPT写悬臂梁挠度计算函数不提任何背景它给我返回了一段基于欧拉-伯努利梁理论的MATLAB风格代码单位默认米制边界条件用的是简支梁还是悬臂梁完全靠猜甚至连输出是挠度还是转角的说明都很模糊。第二个问题我换了一种问法用欧拉-伯努利梁理论计算悬臂梁弯曲。悬臂梁长L2mEI5000 N·m²自由端受集中力P1000N指向-y方向。用有限元法划分20个等长单元固定端节点位移和转角均为零。请给出每个节点的挠度值并将自由端挠度与解析解PL³/(3EI)对比。Python实现使用NumPy。效果天差地别。第二版代码基本一次跑对自由端挠度解析解是0.0533m数值解精确到小数点后四位吻合。差别不在GPT能力而在提问里给了足够多的物理约束。1.2 三段式提示词框架经过大量实验后我现在给GPT布置科学计算编码任务时固定用三段式结构每段都有明确分工。第一段交代物理模型与公式。包括理论名称、几何/材料参数、本构关系、载荷与边界条件。这一段的作用是让GPT在公式层面先和你对齐不至于模型都选错。第二段说明数值实现细节。包括求解域如何离散、单元类型、自由度编号约定、使用什么数值方法处理约束、期望的输入输出格式。第三段要求验证方式与自检用例。这一条最关键必须让GPT在交付代码时同时交付一个独立于主程序的验证脚本用解析解或手算基准检验数值结果。提示这个三段式模板我用了大半年早就背下来了。平时我就是新建一个对话把模板粘贴进去然后改参数和公式描述。强烈建议你也沉淀一份自己的版本而不是每次都临场组织语言。为什么第三段那么重要因为GPT生成代码时有一个很隐蔽的毛病它会极力保持生成内容的内部一致性也就是说主函数里的错误往往会被配套的打印语句、绘图代码掩盖掉——看着输出规律正常实际上物理量已经错了。让它独立写一个和主程序没有耦合的验证函数等于强制它站在对立面去审查自己错误暴露的概率会高很多。2. 实例一2D桁架有限元求解器2.1 任务设定与手算基准这个实例是计算力学入门里最经典的题目之一我用一个最简单的三节点两杆桁架作为基准。节点1坐标(0,0)节点2坐标(2,0)节点3坐标(1,1.5)。杆1连接节点1和3杆2连接节点2和3。两根杆的弹性模量E200GPa截面积A0.01m²。节点1和节点2为固定铰支座节点3受竖直向下的集中力10kN。这个结构足够简单轴力可以用静力平衡手算出来。由于结构对称每根杆分担的竖向力都是5kN。杆与水平方向的夹角θ满足tanθ1.5/11.5sinθ1.5/√(1²1.5²)0.832所以杆件轴力N5/sinθ6.01kN受拉。这个手算结果就是我后面要求GPT验证脚本必须能对上的锚点。2.2 第一轮对话先在公式层面对齐我第一轮提示词不会直接要求它给出完整求解器而是先让它推导并输出单元刚度矩阵的解析表达式确认物理含义一致。提示词大概是这样的我正在用直接刚度法求解一个2D桁架结构。请先写出2D杆单元在整体坐标系下的刚度矩阵表达式用c和s分别表示杆轴线方向余弦和正弦EA为轴向刚度L为杆长。不求代码先列公式。这一步的目的很单纯公式先对齐后面写出来的代码才可能在同一个语义空间里工作。GPT很快给出了标准结果k (EA/L) * [[ c², c·s, -c², -c·s], [ c·s, s², -c·s, -s²], [-c², -c·s, c², c·s], [-c·s, -s², c·s, s²]]这个公式我相信学过有限元的都认识。确认无误后我让它继续写单元刚度矩阵的Python函数只写这一块不写组装。2.3 组装与边界处理的坑接下来提示词进入第二轮和第三轮写全局刚度矩阵组装、处理位移边界条件并求解。我在这里明确了一个自由度编号约定节点i对应全局自由度索引[2i, 2i1]这样节点编号从0开始时组装代码的索引逻辑最直观。GPT给出的组装函数如下import numpy as np def bar_stiffness(x1, y1, x2, y2, E, A): L np.hypot(x2 - x1, y2 - y1) c (x2 - x1) / L s (y2 - y1) / L k E * A / L * np.array([ [c*c, c*s, -c*c, -c*s], [c*s, s*s, -c*s, -s*s], [-c*c, -c*s, c*c, c*s], [-c*s, -s*s, c*s, s*s] ]) return k def assemble_K(n_nodes, nodes, elements, E, A): n_dof 2 * n_nodes K np.zeros((n_dof, n_dof)) for i_global, j_global in elements: xi, yi nodes[i_global] xj, yj nodes[j_global] k bar_stiffness(xi, yi, xj, yj, E, A) idx [2*i_global, 2*i_global1, 2*j_global, 2*j_global1] for a in range(4): for b in range(4): K[idx[a], idx[b]] k[a, b] return K这个阶段GPT表现不错单元刚度矩阵和组装的思路都对。真正容易出问题的是边界条件处理。科学计算里处理固定支座常规有三种方法删行删列法、置大数法、拉格朗日乘子法。我直接要求它用置大数法因为删行删列法在自由度编号后需要重排全局矩阵索引GPT在这种索引重排上出错率很高。置大数法则简单得多在固定自由度对应的对角元上加上一个远大于其它刚度量级的值再把载荷向量对应位置置零for dof in fixed_dofs: K[dof, dof] 1e12 F[dof] 0.0这里有一个很重要的实操细节置大数法加的大数要足够大但又不能大到让矩阵数值条件数恶化。我的经验是取主对角元最大值的1e6到1e8倍而不是随便写一个1e16。曾见过有人用1e16结果双精度浮点下约束自由度上的残余响应反而污染了其他自由度的解。1e12对这个桁架算例是合适的实测求解结果稳定。2.4 让GPT生成对拍脚本从自查变成他查主程序写完后我要求GPT生成验证脚本。提示词是请写一个独立的验证函数不调用主程序用静力平衡条件检查杆件轴力。已知每根杆竖向分力为5kN杆与水平夹角θ满足tanθ1.5检查程序输出的轴力向量是否等于6.01kN受拉为正允许相对误差1e-6。这一步的意义在于GPT如果继续围绕主程序写验证很容易写出主程序输出什么就验证什么的循环论证。要求独立验证后它只能从物理平衡条件出发构造测试这样反而能暴露主程序里的真实错误。第一版验证脚本还真抓到一个问题。GPT在构造单元杆端内力时方向向量是从局部坐标变换回的但它把杆件受拉为正约定弄反了导致验证断言差点失败。我要求它重新检查坐标变换和应力恢复逻辑它才发现局部轴力方向向量组装时少了负号。这种问题单看主程序的位移输出完全看不出来因为节点位移量级和趋势都对但应力恢复环节错了后续基于应力做的设计判断就会全盘出错。2.5 小结一下这个例子的提示词要点这个实例跑下来我提取了三条对后续任务一直有用的经验。第一先对公式再对代码公式阶段GPT出错的概率远低于代码阶段如果公式错了代码阶段纠正的成本高得多。第二自由度编号约定必须显式写进提示词。GPT不会自动知道你心里的编号方案你不说它就用自己训练的默认习惯而这个默认习惯和你的后处理代码很可能不匹配。第三验证函数必须是对立的不能和主程序共享任何内部数据流只能共享物理参数和结果数值。这一点我在多个项目里反复验证过是抓bug最有效的手段。3. 实例二Lennard-Jones势下的原子能量与受力计算3.1 任务设定从势函数到原子受力第二个实例是材料计算里的经典入门任务给定一组原子的坐标用Lennard-Jones势计算体系总势能以及每个原子受到的合力。我用的参数是氩原子的经典值ε0.010323 eVσ3.405 Å。原子数不固定从十几个到上千个都可能有所以代码要考虑效率至少不能写一个完全无脑的纯Python双重循环。这个任务看起来简单实际上坑很多。L-J势的表达式本身很简单但原子受力的推导、力的方向约定、势能截断、截断偏移、pair计数是否重复这些问题一个接一个全是GPT喜欢自由发挥然后出错的地方。3.2 提示词设计先推导再实现再自检我用的提示词完整复述如下请用Lennard-Jones势计算一组原子的总势能和每个原子的受力。L-J势形式为U(r)4ε[(σ/r)^12-(σ/r)^6]ε0.010323 eVσ3.405 Å。坐标系单位为Å能量单位为eV。请先给出力的解析推导式再给出Python实现。实现要求1) 势能函数与受力函数分离2) 考虑截断半径rc2.5σ并在势能中引入U(rc)偏移保证势能在截断处连续3) 力的返回值是逐原子向量并明确说明力的方向约定4) 最后给出一个双原子验证用例。这里每个要求都不是随便写的。先说先推导力的表达式涉及对r求导GPT偶尔会把系数或幂次写错先让它写出解析推导我扫一眼公式就能判断它理解对不对。再说截断偏移这是材料计算的实际需求如果某个体系跑分子动力学模拟势能截断处不连续会产生巨大的伪力所以必须让GPT把U(rc)项算进去。最后说力的方向约定L-J势下两个原子间的作用力是吸引还是排斥方向到底是沿着r_i-r_j还是反方向GPT经常搞混明确要求它说明之后就很难再糊弄过去。GPT给出的受力推导从U(r)4ε[(σ/r)^12-(σ/r)^6]出发对r求导得到dU/dr 24ε/r [2(σ/r)^12 - (σ/r)^6]F_i -∑_j dU/dr_ij · (r_i-r_j)/r_ij等价地写成常见的对拍友好形式F_i 48ε/σ² ∑_j [(σ/r_ij)^14 - 0.5(σ/r_ij)^8] (r_i-r_j)这个形式在文献里很常见优点是把σ/r的主项化成幂次项数值实现时只需要两次除法就能同时得到r^12和r^6的中间量。我让GPT采用这种形式便于后面和实验数据对比。3.3 实现细节截断、偏移与受力更新代码实现部分我要求它必须用NumPy数组操作并明确允许双重循环因为截断半径取2.5σ后每个原子的近邻数在密排体系里也就是十来个双重循环加截断过滤已经完全够用。如果原子数上万再考虑网格分区或近邻列表。GPT生成的代码大致如下import numpy as np def lj_energy_forces(positions, eps, sigma, rc): n_atoms len(positions) U 0.0 F np.zeros_like(positions, dtypefloat) rc2 rc * rc s6_rc (sigma / rc) ** 6 U_shift 4.0 * eps * (s6_rc ** 2 - s6_rc) for i in range(n_atoms): for j in range(i 1, n_atoms): d positions[j] - positions[i] r2 np.dot(d, d) if r2 rc2: continue r np.sqrt(r2) inv_r 1.0 / r s6 (sigma * inv_r) ** 6 s12 s6 * s6 U_pair 4.0 * eps * (s12 - s6) - U_shift U U_pair force_mag 48.0 * eps / (sigma * sigma) * ((sigma * inv_r) ** 14 - 0.5 * (sigma * inv_r) ** 8) force force_mag * d F[i] force F[j] - force return U, F这段代码的标准程度相当高。几个细节处理得很到位一是势能部分减去了U_shift这样在rrc处U_pair恰好为零势能曲线连续二是受力部分用(sigma*inv_r)14和8而不是重新做除法数值上更稳定也更快三是力的更新用了F[i]force和F[j]-force满足牛顿第三定律体系总受力为零。不过这段代码也不是完全没有隐患。有一个非常隐蔽的问题force_mag* d中方向向量d的单位是Åforce_mag的量纲需要仔细检查。实际L-J力公式中48ε/σ²乘以方向向量出来的力单位是eV/Å这个单位在做分子模拟时是可以直接用的。但如果后续要和其他单位制程序对接就需要转换。这种单位一致性问题GPT不会主动提醒必须自己心里有数。3.4 对拍验证双原子手算解析解双原子验证用例是最好用的基准因为L-J势在rσ处的势能恰好为零而力的大小可以手算。两个原子距离rσ时s1U4ε(1-1)0。力的大小F 48ε/σ² [(1)^14 - 0.5(1)^8] × σ方向向量长度 48ε/σ × 0.5 24ε/σ。把ε0.010323 eVσ3.405 Å代进去F24×0.010323/3.4050.0728 eV/Å。方向呢rσ小于平衡距离0.5^(1/6)σ≈1.122σ所以此时两个原子之间是排斥力原子1受到的力指向远离原子2的方向即F_1沿d方向F_2沿-d方向。这个结论可以写进pytest断言里。我要求GPT生成了这样一个验证用例def test_lj_two_atoms_at_sigma(): eps 0.010323 sigma 3.405 rc 2.5 * sigma positions np.array([[0.0, 0.0, 0.0], [sigma, 0.0, 0.0]]) U, F lj_energy_forces(positions, eps, sigma, rc) assert np.isclose(U, 0.0, atol1e-12) expected_force 24.0 * eps / sigma assert np.isclose(np.linalg.norm(F[0]), expected_force, rtol1e-6) assert np.dot(F[0], positions[1] - positions[0]) 0 # 排斥方向 assert np.allclose(F[0] F[1], 0.0, atol1e-12)这个测试本身也是很好的物理锚点训练。GPT在生成主程序时可能注意力全放在实现上忽略方向或系数的错误但一旦要求它写这样的测试它就必须回头审视物理反而能把主程序里潜在的推导错误暴露出来。3.5 材料计算特有的那些看不见的坑L-J实例跑完后我发现材料计算场景下GPT最容易漏的东西非常集中我总结了一个清单截断半径内外的势能连续性处理。不提的话GPT经常不做U_shift修正导致r接近rc时出现伪力尖峰。力的单位与坐标单位的一致性。坐标用Å力自然单位是eV/Å如果程序里混入kJ/mol或者N数值量级会差出十几倍。势能的pair计数问题。我要求的是两个原子间的成对势只算一次所以循环用ij但有些任务里要求的是单个原子的嵌入能或对势总和规则不同必须讲清楚。平衡距离检查L-J势的平衡位置是2^(1/6)σ约1.122σ力在平衡位置为零小于它是排斥大于它是吸引。让GPT写代码前先写出这个区间判断能显著减少方向错误。4. 排错实战GPT代码看起来对、结果错的排查链路4.1 一次真实排错能量曲线系统性偏大的根因有一次我要让GPT生成一个扫描两原子势能曲线的程序横轴是原子间距r从0.8σ到3σ纵轴是总势能U_pair。程序倒是很快写完了但输出曲线和解析解比整体偏大尤其在r接近2σ到2.5σ区间明显不对。我按习惯做三层排查。第一层查单位制输出值解析解都是用eV程序里的eps0.010323也没写错排除。第二层打印每个pair的能量发现r比较大时能量输出比解析值略正r接近rc2.5σ时U_pair应当连续光滑地趋近于0但程序输出在rc附近有一个向正方向的跳变这说明截断处势能没有正确补偿。根因找到了GPT写截断补偿时把U_shift定义成了-U(rc)也就是在pair能量里加了一个正数导致每一个rrc的pair都额外多算了一截正的补偿宏观表现就是总势能系统性偏正。修正方法是把U_shift定义为U(rc)本身也就是一个负数在U_pair里减去它。程序修正后输出曲线和解析解完全重合。这个例子说明一个非常重要的实操经验打印每个pair的贡献比只看总能量高效太多。如果一开始就盯着总能量找bug几个pair错误和正确的可能互相抵消根本定位不到问题。科学计算里一定要养成先局部后整体的排错习惯。4.2 让GPT当代码审查员先列问题清单再改发现截断补偿问题后我没有直接让GPT改代码而是换了一种问法请逐行审查刚才的L-J程序重点检查1) 力与势能的符号关系是否自洽2) 截断处理是否满足rrc处连续3) 能量单位是否一致4) 循环边界是否有重复计数。先不要修改任何代码列出问题清单。这个先列清单再修改的顺序非常关键。GPT如果直接改代码经常会把本来正确的地方也顺手改坏因为它的上下文窗口有限修改过程可能丢失之前的正确设计。让它先列出问题我再挑选真正需要修的点逐一打回去每次只改一处改完立刻跑测试回归。这套流程下来代码的正确性是一次性恢复的不会陷入改A坏B、改B坏C的恶性循环。4.3 最小复现脚本把500行问题化简到20行科学计算代码报错时如果是一个大体系跑出来的错误信息往往淹没在大量无关输出里。我的做法是让GPT生成一个最小复现脚本MRE。提示词是请基于我贴出的完整程序去掉所有文件读写、可视化和日志部分用3个原子、几步迭代构造一个最小脚本触发同样的TypeError错误。为什么是3个原子因为一个原子受力需要至少另一个原子才算有相互作用两个原子能触发力对三个原子则可以暴露周期性、近邻列表和索引边界的问题又足够小到可以手工追踪每个数值。这个办法帮我解决过维度不匹配、近邻列表索引越界、向量化广播错误至少十几次。每次都是MRE一跑错误行号立刻清晰可见。5. 落地工作流把GPT从代码生成器变成推导伙伴5.1 我日常用的提示词模板到这里我想分享一个真实工作中每天都在用的模板写在这里供大家直接抄。背景我在做[材料计算/力学分析]使用Python和NumPy。物理模型模型名称、控制方程、材料参数、几何参数。单位制与数值设定所有量的单位、网格数量/原子数、边界条件、截断参数。任务要求按顺序回答——1) 先列出物理推导步骤和关键公式标注符号含义2) 给出Python实现关键步骤必须注释3) 给出独立验证用例与解析解或手算结果对比。约束不要使用未说明的简化假设如需假设请单独说明并等待确认。这个模板不是万能的但能保证GPT输出结构稳定减少废话最重要的是让物理推导和代码实现分阶段交付便于我快速检查每一步。5.2 多轮对话管理一个任务一个窗口我用GPT做科学计算最大的教训之一就是不要把一个长期项目的所有问题都堆在同一个对话窗口里。GPT的上下文窗口有限越往后它越容易遗忘早期的约定然后开始自由发挥生成的代码里偷偷改变单位、改变边界条件、甚至改变物理模型。我的策略是一个任务一个窗口任务开始前先把该任务相关的所有约定在第一条消息里完整给出。如果中途发现需要补充信息就用注意补充约定...单独成段发出去并且每次让GPT在回复开头复述一遍它当前理解的约定。这样做虽然多花几十个字但能节省后面大量排错时间。5.3 领域知识的最后一公里还是在自己手里说了这么多必须说一点反直觉但很重要的体会GPT写科学计算代码越熟练我越警惕。因为它生成的代码语法几乎永远是对的运行几乎永远不报错但物理可信度必须由人来做最终判断。我自己用的快速检查手段有三层。第一层是量纲分析。看一眼每个公式输出的量纲是否自洽。L-J例子中力的单位是eV/Å还是kJ/mol/nm直接决定后续模拟结果能不能对得上文献。第二层是极限行为。让程序跑几个极端参数r非常大时势能趋近于零、力趋近于零r接近0时势能和力应趋于正无穷平衡距离处力为0。这些极限测试用一行断言就能写却能抓出大量实现错误。第三层是解析解对拍。简单体系总是有解析解或半解析解的哪怕只验证一个最简单的情形也比完全信任程序输出可靠得多。这三层检查配置好后GPT才真正从可能出错的生成器变成了高效的推导与实现伙伴。5.4 最后分享一点个人体会这个系列写到第四篇我对GPT在计算材料与力学编程里的定位越来越清楚。它最大的价值不是替你写代码而是把你从从零写代码里解放出来把精力聚焦到物理判断上。你给它清晰的物理模型、足够的约束边界和可量化的验证手段它就真的能像一位靠谱的结对工程师一样陪你推导公式、实现算法、找出bug。反过来如果你自己心里没有物理基准指望GPT全权负责那最后得到的只会是一堆看起来科学、实际上需要返工的程序。按我的经验每次拿到GPT生成的程序先跑一个解析解对拍再跑一次极限行为检查这两步做到位剩下的问题大多都能在上线前暴露干净。这就是我目前在计算材料写作里的标准流程希望这篇的实例拆解对你有用。

GPT辅助科学计算编程：两个实例拆解提示词设计与验证

相关推荐

P3010 Dividing the Gold 题解：0/1背包转换与方案数DP详解

Pi 1.0 发布：原生MCP与Durable如何重塑终端编程代理

程序员水会生存指南：把无效会议变成高效时间管理

更多资讯

领域特定评估实战：用 Argilla、Distilabel 与 LightEval 构建考试问答评估流水线（smol-course）

Apache Storm 集群安全加固实战：从 OS 层防护到 Kerberos 认证与 ACL 授权

CMake FindOpenCL 模块全解析：从 find_package 到 OpenCL::OpenCL 导入目标

YOLO船舶检测实战：数据集解析与训练避坑指南

题解：洛谷 P14361 [CSP-S 2025] 社团招新

U-Boot Kbuild深度解析：从零构建RV1106移植的四大核心步骤

今日更新

PHP htmlentities()函数用法讲解

AI Agent工程实战：从七要素到七个决策点的系统设计指南

AI入局芯片设计验证：验证收敛50倍加速的真相与应对

本周精选

CMOS设计三大物理效应：体效应、沟道长度调制与亚阈值导电

微信小程序在线阅读毕设全解析：Java后端+MySQL+避坑实战

Android音频问题五层分析法：从App到Codec的系统化诊断

本月热门

我发现了一个新思路：用 Remotion + Claude Code 像写代码一样自动化生成短视频

Windows下 Codex 中 Chrome 和 Computer Use 插件不可用问题排查及解决参考方式：TaoToken 统一 Key 配置与验证

2026 大模型集体涨价：用 Python 做企业 Token 成本测算与选型避坑（附配置）

读完文章，想聊聊您的网站？