lambda算子 1.b

Lambda演算基础
本文深入探讨了Lambda演算中的自由与有界标识符概念,解释了标识符如何在不同的上下文中区分,以及如何通过Alpha转换和Beta简化进行计算。此外还介绍了如何避免变量冲突等问题。
上上周就快写完这篇时,IE突然当掉,写的东西烟消云散。俺也元气大伤。这次吸取教训,不用CSDN的在线工具写了。这样的坏处是把文章拷贝到CSDN时,格式难免出错,还得手工调整一下。CSDN什么时候可以实现confluence的在线备份功能呢?还是继续八卦lambda:

自由 vs 有界标识符

标识符和变量其实是一个意思。我记得国内教材里很少用标识符这个说法。不过既然原作者用这个说法,我就跟着用了。上次说到Currying解决了如何处理多参数的问题。在讨论怎么使用lambda前,我们还要解决一个细微但重要的语法问题:封闭(closure),或者叫完全有界(complete bounding)。这里的有界和一阶谓词逻辑里的有界没有本质区别,对一阶谓词逻辑熟悉的老大们可以放心跳过。其实有界涉及的定义很直观,我们看一个例子先。假设我们有一个函数lambda x y. (lambda y. y + 2) + x + y +z,lambda y. y+2里的y和它后面的y是不是一样的呢?显然它们是不一样的。为了处理这种区别,我们引入了有界。当一个lambda表达式被计算时,它不能处理无界的标识符。当一个标识符出现在一个lambda表达式的参数里,并且被包含在这个lambda表达式里,我们就可以说这个标识符有界。如果一个标识符没有被包含在任何一个表达式里,我们就叫它为自由变量。比如说,上面那个lambda表达式里,x 出现在lambda x y .(....)里,所以它是有界的变量,它的包含环境(enclosing context,用“语境”或者“上下文”怎么听怎么别扭,好像俺是《读书》的御用作者似的。:-D)是整个lambda表达式。lambda y. y+2里的y也是有界的,但它的包含环境是lambda y. y+2。标识符z没有出现在包含它的表达式的参数列表里,所以是自由变量。再举几个例子:

  1. lambda x . plus x y: 这个表达式里,"y"和"plus"都是自由变量,因为它们不是任何包含它们的表达式的参数。x有界,因为它被包含在plus x y里,而plus x y的参数有x。
  2. lambda x y.y x: 这个表达式里,x和y都有界,因为它们是这个表达式的参数。
  3. lambda y . (lambda x . plus x y): 在内嵌的表达式lambda x. plus x y里,y和 plus 是自由变量而x是有界变量。在整个lambda表达式里,x和y都有界:x在内嵌表达式界内,而y在整个表达式界内。plus仍然自由。

我们用"free(x)"来代表表达式x里所有自由变量的集合。

一个lambda表达式完全合法仅当它的所有变量都有界。不过当我们考查某个复杂表达式里的子表达式且不考虑上下文时,那些子表达式可以有自由变量-其实确保正确处理那些子表达式里的自由变量非常重要。

 

Lambda 算子计算规则

其实真正的规则就俩:alpha和beta。Alpha规则又叫转换(conversion)规则,而beta规则又叫简化(reduction)规则。

Alpha转换

这个充满了《星际迷航》味道的规则其实就是重命名操作。它无非是说变量名不重要:给定任何一个lambda表达式,我们可以任意改变参数的名字,只要我们相应地改变这些对应这些参数的变量名字。

比如说,我们有如下表达式:
lambda x . if (= x 0) then 1 else x^2

我们通过alpha规则把X改成Y(写作alpha[x/y], 和逻辑里的变量替换一个写法),于是得到:
lambda y . if (= y 0) then 1 else y^2

Alpha操作完全不影响lambda表达式的意义。不过我们后面会发现,这个操作很重要,因为它让我们能够实现诸如递归的操作。

Beta简化

Beta简化就有意思了。我们只需要这一个规则,就可以让lamdba算子实现一台计算机能做的任何计算。透过纷繁的表象,我们会发现事情的本质往往出人意料地清晰而简单。删繁为简,恰是数学魅力所在。

Beta规则无非是说,应用一个函数(也就是lambda表达式。一个意思)等价于通过把函数体内有界的变量替换成应用里对应参数的实际值来替换原来的函数。听上去有些拗口(呵呵,其实原文更拗口),但当你看一个例子就知道它其实很简单:

假设我们要应用一个函数:"(lambda x . x + 1) 3"。Beta规则说,我们可以替换整个表达式,把函数体(也就是“x+1”)里的参数对应的x替换成实际的值3。所以最后的结果是“3+1”。

再来一个稍微复杂点的例子:
lambda y . (lambda x . x + y)) q

这个表达式有意思,因为应用了这个表达式后,我们可以得到另外一个表达式。也就是说,它是一个生成表达式的表达式(说到这里,玩儿动态语言的老大们可以笑了,玩儿C/C++/Java的老大们可以流口水了)。当我们对这个表达式应用Beta简化时,我们把所有对应参数y的变量替换成实际的值q。所以结果是"lambda x, x+q"。

再来一个例子:
"(lambda x y. x y) (lambda z . z * z) 3". 这个带两个参数的函数把第一个参数应用到第二个参数上。当我们计算它的值时,我们把第一个lambda表达式里的变量x换成lambda z. z * z, 再把变量y换成3,得到(lambda z. z * z) 3。对该结果应用Beta简化,我们得到3 * 3。

Beta的严格定义如下:

lambda x . B e = B[x := e] if free(e) /subset free(B[x := e]

这个定义末尾的条件,"if free(e) /subset free(B[x:=e])"道出了我们需要Alpha转换的原因:仅当beta化简不会引起有界变量和自由变量的冲突时,我们可以实施Beta化简。如果一个变量“z”是"e"里的自由变量,那我们得保证beta化简不会让"z"变成有界变量。如果B里的有界变量和”e"里的自由变量重名,我们必须先用Alpha转换,是的重名的变量不再重名。形式化定义不够直观,直观描述又不够简洁。还是来个例子漱漱口:

给定一个表达式,lambda z. lambda x. x+z. 假设我们要应用这个表达式:
(lambda z . (lambda x . x + z)) (x + 2)

在实际参数"(x + 2)"里,x是自由变量。但x不是表达式lambda x. x+z的自由变量。也就是说,free(e) /subset free(B[ x:=e])不成立。如果我们打破Beta简化的规则,直接开始Beta简化,便会得到: lambda x . x + x + 2

"x+2"里自由变量,x,现在变得有界了。如果我们把结果应用到3上:(lambda x. x+2+2) 3,我们得到3 + 3 + 2。

如果我们按正常程序办事呢?

应用 alpha[x/y]: lambda z . (lambda y . y+z)) (x + 2)
应用 beta: lambda y . y + x + 2) 3
再次应用beta: 3 + x + 2.

"3+x+2" 和 "3+3+2" 很不一样哈!

规则就这些了。我们还可以选择性地加一个所谓的Eta-化简,不过它不是必需的。我们就此跳过。我们讨论的这套系统已经是图灵完备的计算体系。那这套系统到底有什么用嗫?到底怎样才能让这套系统变得真正有用嗫?嗯,要说明这些问题,我们得先定义一些基本的函数,以便我们做算术,条件测试,递归,等等。这些会在以后的帖子里谈到。

我们也还没有谈到适合lambda算子的模型(Good Math Bad Math的作者在这里这里讨论了模型)。模型也是很重要的东西。逻辑学家用了好几年时间研究lambda算子,才搞出一个完备的模型。而且早先时候,尽管lambda算子看起来没错,为它制订模型的工作却失败了。这在当时极为引人关注。要知道,毕竟一个系统没有有效的模型就没有实际的意义。

查找错误 function [psinew, phinew] = new2D1(h, dt) % 二维非线性薛定谔方程求解器 L = 10; a = 1; % 空间域半长,波速 T = 1; % 终止时间 M = round(2*L/h); % 每维空间网格数 x = (-L:h:L)'; x = x(1:end-1); % 周期边界处理 [x, y] = meshgrid(x, x); % 二维网格 % 非线性参数 a1 = 0.5; b1 = -2; c1 = -1; a2 = 0.5; b2 = -3; c2 = -1; % 2-stage Gauss系数 c1g = 0.5 - sqrt(3)/6; c2g = 0.5 + sqrt(3)/6; a11 = 1/4; a12 = 1/4 - sqrt(3)/6; a21 = 1/4 + sqrt(3)/6; a22 = 1/4; b1g = 0.5; b2g = 0.5; %% 计算二维拉普拉斯算子的特征值 kk = 0:M-1; theta = 2*pi*kk/M; lambda_A = (150 + 32*cos(theta) - 2*cos(2*theta)) / 180; coef_S = 1/(60*h^2); lambda_S = (-126 + 128*cos(theta) - 2*cos(2*theta)) * coef_S; lambda_L = lambda_S ./ lambda_A; % 二维特征值矩阵 [KX, KY] = meshgrid(lambda_L, lambda_L); lambda_B = KX + KY; % 二维拉普拉斯算子特征值 %% 预计算FFT求解器 invD = zeros(4, 4, M, M); % 每个频率点的4x4矩阵 for i = 1:M for j = 1:M lam = lambda_B(i, j); D_k = [1, 0, -a^2*dt*a11*lam, -a^2*dt*a12*lam; 0, 1, -a^2*dt*a21*lam, -a^2*dt*a22*lam; -dt*a11, -dt*a12, 1, 0; -dt*a21, -dt*a22, 0, 1]; invD(:, :, i, j) = inv(D_k); % 预计算逆矩阵 end end %% 初始条件 - 二维高斯波包 r = sqrt(x.^2 + y.^2); exact_psi = @(x, y, t) 4*atan(exp(3 - 5*r)); exact_phi = @(x, y, t) atan(exp(3 - 5*r)); psi = exact_psi(x, y, 0); psi0 = psi; phi = exact_phi(x, y, 0); phi0 = phi; u = zeros(size(psi)); v = zeros(size(phi)); %% 初始质量 M0 M0 = h^2 * sum(sum(abs(psi0).^2 + abs(phi0).^2)); % 初始能量 E0 term1 = h^2 * sum(sum(u.^2)); term2 = -a^2 * h^2 * sum(sum(real(ifft2(lambda_B .* fft2(psi0))) .* psi0)); term3 = a1 * h^2 * sum(sum(psi0.^2)); term4 = h^2 * sum(sum(v.^2)); term5 = -a^2 * h^2 * sum(sum(real(ifft2(lambda_B .* fft2(phi0))) .* phi0)); term6 = a2 * h^2 * sum(sum(phi0.^2)); term7 = (c2*b1/(4*c1)) * h^2 * sum(sum(psi0.^4)); term8 = (b2/4) * h^2 * sum(sum(phi0.^4)); term9 = (c2/2) * h^2 * sum(sum((psi0.*phi0).^2)); E0 = (1/2) * ((c2/c1)*(term1 + term2 + term3) + term4 + term5 + term6) + term7 + term8 + term9; Nt = round(T/dt); % 初始化守恒量记录 M_arr = zeros(Nt+1, 1); E_arr = zeros(Nt+1, 1); abs_err_M = zeros(Nt+1, 1); abs_err_E = zeros(Nt+1, 1); M_arr(1) = M0; E_arr(1) = E0; abs_err_M(1) = 0; abs_err_E(1) = 0; %% 初始化 prev_u = u; prev_v = v; prev_psi = psi; prev_phi = phi; prev_psi1 = psi; prev_psi2 = psi; prev_phi1 = phi; prev_phi2 = phi; prev_u1 = u; prev_u2 = u; prev_v1 = v; prev_v2 = v; prev_t = 0; %% 时间层循环 for n = 1:Nt %% 使用预测龙格库塔求高阶预估值 current_t = n*dt; u = prev_u; u1 = prev_u1; u2 = prev_u2; v = prev_v; v1 = prev_v1; v2 = prev_v2; psi = prev_psi; psi1 = prev_psi1; psi2 = prev_psi2; phi = prev_phi; phi1 = prev_phi1; phi2 = prev_phi2; % 使用前一时间步的值进行三点插值 t_nodes = [prev_t, prev_t + c1g*dt, prev_t + c2g*dt]; U_nodes = [u, u1, u2]; V_nodes = [v, v1, v2]; psi_nodes = [psi, psi1, psi2]; phi_nodes = [phi, phi1, phi2]; % 当前阶段点时间 t_stage1 = current_t + c1g*dt; t_stage2 = current_t + c2g*dt; % 内联拉格朗日插值计算 u1star0 = lagrange3_2D(t_stage1, t_nodes, u, u1, u2); u2star0 = lagrange3_2D(t_stage2, t_nodes, u, u1, u2); v1star0 = lagrange3_2D(t_stage1, t_nodes, v, v1, v2); v2star0 = lagrange3_2D(t_stage2, t_nodes, v, v1, v2); psi1star0 = lagrange3_2D(t_stage1, t_nodes, psi, psi1, psi2); psi2star0 = lagrange3_2D(t_stage2, t_nodes, psi, psi1, psi2); phi1star0 = lagrange3_2D(t_stage1, t_nodes, phi, phi1, phi2); phi2star0 = lagrange3_2D(t_stage2, t_nodes, phi, phi1, phi2); % 迭代初值 u1n_s = u1star0; psi1n_s = psi1star0; u2n_s = u2star0; psi2n_s = psi2star0; v1n_s = v1star0; phi1n_s = phi1star0; v2n_s = v2star0; phi2n_s = phi2star0; iterw = 0; err = inf; % 迭代循环 while iterw < 50 && err > 1e-14 % 非线性项 p1n_s = a1*psi1n_s + b1*psi1n_s.^3 + c1*psi1n_s.*phi1n_s.^2; p2n_s = a1*psi2n_s + b1*psi2n_s.^3 + c1*psi2n_s.*phi2n_s.^2; q1n_s = a2*phi1n_s + b2*phi1n_s.^3 + c2*phi1n_s.*psi1n_s.^2; q2n_s = a2*phi2n_s + b2*phi2n_s.^3 + c2*phi2n_s.*psi2n_s.^2; % 列方程组求解 bw = cat(3,u - dt*a11*p1n_s - dt*a12*p2n_s,... u - dt*a21*p1n_s - dt*a22*p2n_s,... psi,... psi); bo = cat(3, v - dt*a11*q1n_s - dt*a12*q2n_s, ... v - dt*a21*q1n_s - dt*a22*q2n_s, ... phi, ... phi); % 使用FFT求解器 xw = solve_fft2(bw, invD, M); xo = solve_fft2(bo, invD, M); % 更新迭代值 u1n_new = xw(:,:,1); u2n_new = xw(:,:,2); psi1n_new = xw(:,:,3); psi2n_new = xw(:,:,4); v1n_new = xo(:,:,1); v2n_new = xo(:,:,2); phi1n_new = xo(:,:,3); phi2n_new = xo(:,:,4); % 计算误差 err_u1n = max(abs(u1n_new(:) - u1n_s(:))); err_u2n = max(abs(u2n_new(:) - u2n_s(:))); err_psi1n = max(abs(psi1n_new(:) - psi1n_s(:))); err_psi2n = max(abs(psi2n_new(:) - psi2n_s(:))); err_v1n = max(abs(v1n_new(:) - v1n_s(:))); err_v2n = max(abs(v2n_new(:) - v2n_s(:))); err_phi1n = max(abs(phi1n_new(:) - phi1n_s(:))); err_phi2n = max(abs(phi2n_new(:) - phi2n_s(:))); err = max([err_u1n, err_u2n, err_psi1n, err_psi2n,... err_v1n, err_v2n, err_phi1n, err_phi2n]); % 更新迭代值 u1n_s = u1n_new; u2n_s = u2n_new; psi1n_s = psi1n_new; psi2n_s = psi2n_new; v1n_s = v1n_new; v2n_s = v2n_new; phi1n_s = phi1n_new; phi2n_s = phi2n_new; iterw = iterw + 1; end % 输出最终迭代值 psi1_star = psi1n_s; psi2_star = psi2n_s; phi1_star = phi1n_s; phi2_star = phi2n_s; %% 计算非线性项 p1n_star = a1*psi1_star + b1*psi1_star.^3 + c1*psi1_star.*phi1_star.^2; p2n_star = a1*psi2_star + b1*psi2_star.^3 + c1*psi2_star.*phi2_star.^2; q1n_star = a2*phi1_star + b2*phi1_star.^3 + c2*phi1_star.*psi1_star.^2; q2n_star = a2*phi2_star + b2*phi2_star.^3 + c2*phi2_star.*psi2_star.^2; % 右端项 bb0 = cat(3, a^2*real(ifft2(lambda_B.*fft2(psi))), ... a^2*real(ifft2(lambda_B.*fft2(psi))), ... u, u); bb1 = cat(3, -p1n_star, zeros(M), zeros(M), zeros(M)); bb2 = cat(3, zeros(M), -p2n_star, zeros(M), zeros(M)); bb3 = cat(3, -psi1_star, zeros(M), zeros(M), zeros(M)); bb4 = cat(3, zeros(M), -psi2_star, zeros(M), zeros(M)); % 使用FFT求解器 x0 = solve_fft2(bb0, invD, M); x1 = solve_fft2(bb1, invD, M); x2 = solve_fft2(bb2, invD, M); x3 = solve_fft2(bb3, invD, M); x4 = solve_fft2(bb4, invD, M); % 右端项 yb0 = cat(3, a^2*real(ifft2(lambda_B.*fft2(phi))), ... a^2*real(ifft2(lambda_B.*fft2(phi))), ... v, v); yb1 = cat(3, -q1n_star, zeros(M), zeros(M), zeros(M)); yb2 = cat(3, zeros(M), -q2n_star, zeros(M), zeros(M)); yb3 = cat(3, -phi1_star, zeros(M), zeros(M), zeros(M)); yb4 = cat(3, zeros(M), -phi2_star, zeros(M), zeros(M)); % 使用FFT求解器 y0 = solve_fft2(yb0, invD, M); y1 = solve_fft2(yb1, invD, M); y2 = solve_fft2(yb2, invD, M); y3 = solve_fft2(yb3, invD, M); y4 = solve_fft2(yb4, invD, M); % 提取各部分 k1_1_0 = x0(:,:,1); k1_1_lamda1 = x1(:,:,1); k1_1_lamda2 = x2(:,:,1); k1_1_eita1 = x3(:,:,1); k1_1_eita2 = x4(:,:,1); k1_2_0 = x0(:,:,2); k1_2_lamda1 = x1(:,:,2); k1_2_lamda2 = x2(:,:,2); k1_2_eita1 = x3(:,:,2); k1_2_eita2 = x4(:,:,2); k2_1_0 = x0(:,:,3); k2_1_lamda1 = x1(:,:,3); k2_1_lamda2 = x2(:,:,3); k2_1_eita1 = x3(:,:,3); k2_1_eita2 = x4(:,:,3); k2_2_0 = x0(:,:,4); k2_2_lamda1 = x1(:,:,4); k2_2_lamda2 = x2(:,:,4); k2_2_eita1 = x3(:,:,4); k2_2_eita2 = x4(:,:,4); k3_1_0 = y0(:,:,1); k3_1_lamda1 = y1(:,:,1); k3_1_lamda2 = y2(:,:,1); k3_1_eita1 = y3(:,:,1); k3_1_eita2 = y4(:,:,1); k3_2_0 = y0(:,:,2); k3_2_lamda1 = y1(:,:,2); k3_2_lamda2 = y2(:,:,2); k3_2_eita1 = y3(:,:,2); k3_2_eita2 = y4(:,:,2); k4_1_0 = y0(:,:,3); k4_1_lamda1 = y1(:,:,3); k4_1_lamda2 = y2(:,:,3); k4_1_eita1 = y3(:,:,3); k4_1_eita2 = y4(:,:,3); k4_2_0 = y0(:,:,4); k4_2_lamda1 = y1(:,:,4); k4_2_lamda2 = y2(:,:,4); k4_2_eita1 = y3(:,:,4); k4_2_eita2 = y4(:,:,4); %% 求解优化问题(简化版) % 在实际应用中,这里需要求解一个非线性系统来获得最优参数 % 为简化,我们使用固定值 lamda1 = 1; lamda2 = 1; eita1 = 0; eita2 = 0; kesi1 = 1; kesi2 = 1; %% 组装得到K1(1,2),K2(1,2),K3(1,2),K4(1,2) k1_1 = k1_1_0 + lamda1.*k1_1_lamda1 + lamda2.*k1_1_lamda2 + eita1.*k1_1_eita1 + eita2.*k1_1_eita2; k1_2 = k1_2_0 + lamda1.*k1_2_lamda1 + lamda2.*k1_2_lamda2 + eita1.*k1_2_eita1 + eita2.*k1_2_eita2; k2_1 = k2_1_0 + lamda1.*k2_1_lamda1 + lamda2.*k2_1_lamda2 + eita1.*k2_1_eita1 + eita2.*k2_1_eita2; k2_2 = k2_2_0 + lamda1.*k2_2_lamda1 + lamda2.*k2_2_lamda2 + eita1.*k2_2_eita1 + eita2.*k2_2_eita2; k3_1 = k3_1_0 + lamda1.*k3_1_lamda1 + lamda2.*k3_1_lamda2 + eita1.*k3_1_eita1 + eita2.*k3_1_eita2; k3_2 = k3_2_0 + lamda1.*k3_2_lamda1 + lamda2.*k3_2_lamda2 + eita1.*k3_2_eita1 + eita2.*k3_2_eita2; k4_1 = k4_1_0 + lamda1.*k4_1_lamda1 + lamda2.*k4_1_lamda2 + eita1.*k4_1_eita1 + eita2.*k4_1_eita2; k4_2 = k4_2_0 + lamda1.*k4_2_lamda1 + lamda2.*k4_2_lamda2 + eita1.*k4_2_eita1 + eita2.*k4_2_eita2; %% n+1层的解 unew = u + dt*b1g*k1_1 + dt*b2g*k1_2; psinew = psi + dt*b1g*k2_1 + dt*b2g*k2_2; vnew = v + dt*b1g*k3_1 + dt*b2g*k3_2; phinew = phi + dt*b1g*k4_1 + dt*b2g*k4_2; % 计算当前时间步的质量和能量 M_current = h^2 * sum(sum(abs(psinew).^2 + abs(phinew).^2)); term1 = h^2 * sum(sum(vnew.^2)); term2 = -a^2 * h^2 * sum(sum(real(ifft2(lambda_B .* fft2(psinew))) .* psinew)); term3 = a1 * h^2 * sum(sum(psinew.^2)); term4 = h^2 * sum(sum(unew.^2)); term5 = -a^2 * h^2 * sum(sum(real(ifft2(lambda_B .* fft2(phinew))) .* phinew)); term6 = a2 * h^2 * sum(sum(phinew.^2)); term7 = (c2*b1/(4*c1)) * h^2 * sum(sum(psinew.^4)); term8 = (b2/4) * h^2 * sum(sum(phinew.^4)); term9 = (c2/2) * h^2 * sum(sum((psinew.*phinew).^2)); E_current = (1/2)*((c2/c1)*(term1 + term2 + term3) + term4 + term5 + term6) + term7 + term8 + term9; % 记录守恒量及误差 M_arr(n+1) = M_current; E_arr(n+1) = E_current; abs_err_M(n+1) = abs(M_current - M0); abs_err_E(n+1) = abs(E_current - E0); % 更新历史变量 prev_psi1 = psi + dt*(a11*k2_1 + a12*k2_2); prev_psi2 = psi + dt*(a21*k2_1 + a22*k2_2); prev_phi1 = phi + dt*(a11*k4_1 + a12*k4_2); prev_phi2 = phi + dt*(a21*k4_1 + a22*k4_2); prev_u1 = u + dt*(a11*k1_1 + a12*k1_2); prev_u2 = u + dt*(a21*k1_1 + a22*k1_2); prev_v1 = v + dt*(a11*k3_1 + a12*k3_2); prev_v2 = v + dt*(a21*k3_1 + a22*k3_2); % 更新n+1层的值 prev_u = unew; prev_psi = psinew; prev_v = vnew; prev_phi = phinew; prev_t = n*dt; % 每10步显示进度 if mod(n, 10) == 0 fprintf('完成时间步:%d/%d, 质量误差:%.4e, 能量误差:%.4e\n', ... n, Nt, abs_err_M(n+1), abs_err_E(n+1)); end end %% 结果可视化 % 绘制守恒量误差图 figure; semilogy((0:Nt)*dt, abs_err_M, 'b', 'LineWidth', 1.5); xlabel('t'); ylabel('绝对误差 (log scale)'); title('质量守恒误差'); grid on; figure; semilogy((0:Nt)*dt, abs_err_E, 'r', 'LineWidth', 1.5); xlabel('t'); ylabel('绝对误差 (log scale)'); title('能量守恒误差'); grid on; % 绘制初始和最终波形 figure; subplot(2,2,1); surf(x, y, psi0, 'EdgeColor', 'none'); title('初始 \psi'); xlabel('x'); ylabel('y'); subplot(2,2,2); surf(x, y, psinew, 'EdgeColor', 'none'); title('最终 \psi'); xlabel('x'); ylabel('y'); subplot(2,2,3); surf(x, y, phi0, 'EdgeColor', 'none'); title('初始 \phi'); xlabel('x'); ylabel('y'); subplot(2,2,4); surf(x, y, phinew, 'EdgeColor', 'none'); title('最终 \phi'); xlabel('x'); ylabel('y'); end %% 二维FFT求解器函数 function x = solve_fft2(b, invD, N) % 输入b是N×N×4的三维数组 % 输出x是N×N×4的三维数组 % FFT变换 B1 = fft2(b(:,:,1)); B2 = fft2(b(:,:,2)); B3 = fft2(b(:,:,3)); B4 = fft2(b(:,:,4)); % 预分配解的空间 X1 = zeros(N, N); X2 = zeros(N, N); X3 = zeros(N, N); X4 = zeros(N, N); % 频域求解 for i = 1:N for j = 1:N rhs = [B1(i,j); B2(i,j); B3(i,j); B4(i,j)]; sol = invD(:,:,i,j) * rhs; X1(i,j) = sol(1); X2(i,j) = sol(2); X3(i,j) = sol(3); X4(i,j) = sol(4); end end % IFFT变换回物理空间 x1 = real(ifft2(X1)); x2 = real(ifft2(X2)); x3 = real(ifft2(X3)); x4 = real(ifft2(X4)); % 组合解 x = cat(3, x1, x2, x3, x4); end %% 二维三点拉格朗日插值 function y = lagrange3_2D(t, t_nodes, y0, y1, y2) t0 = t_nodes(1); t1 = t_nodes(2); t2 = t_nodes(3); L0 = ((t - t1).*(t - t2)) ./ ((t0 - t1).*(t0 - t2)); L1 = ((t - t0).*(t - t2)) ./ ((t1 - t0).*(t1 - t2)); L2 = ((t - t0).*(t - t1)) ./ ((t2 - t0).*(t2 - t1)); y = y0.*L0 + y1.*L1 + y2.*L2; end
最新发布
09-07
该代码实现了一个二维非线性薛定谔方程的求解器,但存在一些潜在的问题和错误,以下是详细的分析和建议: ### 1. **空间网格生成问题** ```matlab x = (-L:h:L)'; x = x(1:end-1); % 周期边界处理 [x, y] = meshgrid(x, x); % 二维网格 ``` - **问题**:`x = (-L:h:L)'` 生成的点数可能不是整数,导致网格点不均匀。 - **建议**:使用 `linspace` 或确保 `L/h` 是整数,并通过 `M = round(2*L/h)` 保证整数个点。 ### 2. **特征值计算错误** ```matlab lambda_A = (150 + 32*cos(theta) - 2*cos(2*theta)) / 180; coef_S = 1/(60*h^2); lambda_S = (-126 + 128*cos(theta) - 2*cos(2*theta)) * coef_S; lambda_L = lambda_S ./ lambda_A; ``` - **问题**:`lambda_A` 和 `lambda_S` 的公式可能不正确,尤其是用于二维拉普拉斯算子的近似。这些公式可能没有正确反映空间离散化后的特征值。 - **建议**:确认这些公式是否符合所使用的数值离散方法(如有限差分或谱方法)。 ### 3. **能量计算错误** ```matlab term2 = -a^2 * h^2 * sum(sum(real(ifft2(lambda_B .* fft2(psi0))) .* psi0)); ``` - **问题**:`ifft2(lambda_B .* fft2(psi0))` 计算的是拉普拉斯作用后的结果,但 `real` 函数可能不必要,且 `psi0` 的点乘可能不正确。 - **建议**:直接使用 `fft2` 和 `ifft2` 计算动能项,避免手动处理实部。 ### 4. **迭代求解器收敛性问题** ```matlab while iterw < 50 && err > 1e-14 ... err = max([err_u1n, err_u2n, err_psi1n, err_psi2n,... err_v1n, err_v2n, err_phi1n, err_phi2n]); end ``` - **问题**:迭代次数上限为 50,但收敛条件 `err > 1e-14` 可能过于严格,导致无法收敛。 - **建议**:调整收敛阈值或增加迭代次数上限,或者添加不收敛时的警告提示。 ### 5. **非线性项计算错误** ```matlab p1n_s = a1*psi1n_s + b1*psi1n_s.^3 + c1*psi1n_s.*phi1n_s.^2; ``` - **问题**:非线性项的表达式可能与物理模型不符,尤其是 `psi1n_s.*phi1n_s.^2` 这一项。 - **建议**:检查非线性项是否与物理模型一致,确保数学表达式正确。 ### 6. **时间步长与稳定性** ```matlab Nt = round(T/dt); ``` - **问题**:`dt` 可能过大,导致数值不稳定。 - **建议**:确保 `dt` 满足 CFL 条件,或者在代码中添加自动调整时间步长的功能。 ### 7. **函数 `lagrange3_2D` 插值问题** ```matlab function y = lagrange3_2D(t, t_nodes, y0, y1, y2) ``` - **问题**:该函数是否正确处理了二维数据? - **建议**:检查 `lagrange3_2D` 是否适用于矩阵输入,尤其是 `y0`, `y1`, `y2` 为矩阵时的运算是否正确。 ### 8. **输出变量未定义** ```matlab function [psinew, phinew] = new2D1(h, dt) ``` - **问题**:如果 `psinew` 和 `phinew` 在某些情况下未被正确赋值,可能导致输出错误。 - **建议**:确保所有可能的执行路径都为 `psinew` 和 `phinew` 赋值。 ### 9. **可视化部分问题** ```matlab surf(x, y, psi0, 'EdgeColor', 'none'); ``` - **问题**:`x` 和 `y` 的维数是否与 `psi0` 匹配? - **建议**:确保 `x` 和 `y` 的大小与 `psi0` 一致,否则可能导致绘图错误。 ### 10. **变量命名不一致** - **问题**:变量名如 `psi`, `psinew`, `psi0` 等容易混淆。 - **建议**:统一命名规则,增加注释以提高可读性。 --- ###
评论 6
成就一亿技术人!
拼手气红包6.0元
还能输入1000个字符
 
红包 添加红包
表情包 插入表情
 条评论被折叠 查看
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值