电磁场有限元方法EX2.2-里兹法求解泊松方程控制的边值问题

电磁场有限元方法EX2.2-里兹法求解泊松方程控制的边值问题

简单学习一下有限元法的基础理论，书本为电磁场有限元经典教材：
THE FINITE ELEMENT METHOD IN ELECTROMAGNETICS, JIAN-MING JIN

1、问题定义

1.1、英文描述

1.2、中文描述

给定由波动方程定义的边值问题：

$$\frac{{\mathrm{d}}^2\varphi }{\mathrm{d}{x}^2}+2\varphi =0(0<x<1)$$

边界条件：

$${\varphi |}_{x=0}=0\text{和}{\varphi |}_{x=1}=1$$

证明该问题的泛函$F$为：

$$F(\tilde{\varphi })=\frac{1}{2}{\int }_0^1\left[{\left(\frac{\mathrm{d}\tilde{\varphi }}{\mathrm{d}x}\right)}^2-2{\tilde{\varphi }}^2\right]\mathrm{d}x$$

应用里茨方法，使用试探函数：

$$\tilde{\varphi }(x)=x+\sum _{n=1}^N{c}_n\sin (n\pi x)(\text{取 }N=2)$$

求近似解，并将结果与精确解：

$$\varphi =\frac{\sin (\sqrt{2}x)}{\sin \sqrt{2}}$$

进行比较。

2、解题

2.1、解题思路

这题比较简单，使用经典里兹方法，将试探函数$\tilde{\varphi }(x)$带入泛函$F$中，然后对泛函$F$求c1和c2的偏导即可（$F$对c1，c2的偏导为0，以此求得极小值）。变分公式的证明见最后。

下面使用matlab实现代码思路

2.2、Matlab代码实现

2.2.1 定义试探函数和变分问题

定义试探函数和变分问题，并将试探函数带入变分：

%% Solve
syms x c1 c2
% define the trail function and variational functional
phi = x+c1*sin(pi*x) +c2*sin(2*pi*x);
F=0.5*int((diff(phi, x))^2-2*(phi^2),x,0,1);
disp('F=');
pretty(simplify(F))

2.2.2 变分F对c1、c2变量求偏导

% Find the extreme value of the functional. calculate the partial derivative of c
diff_c1 = diff(F, c1);  % ∂F/∂c1
diff_c2 = diff(F, c2);  % ∂F/∂c2

2.2.3 F对c1、c2变量偏导为0，建立方程并求解

% Solve the system of equations for variables [c1, c2]
equations = [diff_c1 == 0, diff_c2 == 0];
solutions = solve(equations, [c1, c2]);
disp('Solution:');
% Extract all solutions for c1 and c2
c1_val = double(solutions.c1);
c2_val = double(solutions.c2);
% Display results (both symbolic and numeric forms)
disp(['c1 = ', num2str(c1_val)]);
disp(['c2 = ', num2str(c2_val)]);

2.2.4 对比理论解和数值解，并绘图

phi_curve=subs(phi, {c1, c2}, {c1_val, c2_val});
phi_curve_exact=sin(sqrt(2)*x)/sin(sqrt(2));%% Draw figure
figure(1)
fplot(phi_curve, [0, 1], ...'Color', [0, 0.4470, 0.7410], ... % MATLAB经典蓝'LineStyle', '--', ...             % 虚线'LineWidth', 2.5, ...'Marker', 'o', ...                 % 圆形标记'MarkerSize', 7, ...'MarkerFaceColor', 'auto', ...     % 自动填充'MarkerEdgeColor', 'k', ...        % 黑色边框'DisplayName', sprintf('Numerical: $\\phi(x) = x + %.2f\\sin(\\pi x) + %.2f\\sin(2\\pi x)$', c1_val, c2_val));hold on
fplot(phi_curve_exact, [0, 1], ...'Color', [0.8500, 0.3250, 0.0980], ... % MATLAB经典红'LineStyle', '-', ...                   % 实线'LineWidth', 2.5, ...'Marker', 's', ...                      % 正方形标记'MarkerSize', 8, ...'MarkerFaceColor', 'auto', ...           % 自动填充'DisplayName', 'Exact: $\frac{\sin(\sqrt{2} x)}{\sin(\sqrt{2})}$');xlabel('Position (x)', 'FontSize', 14);
ylabel('Function Value \phi(x)', 'FontSize', 14);
title('Comparison of Numerical and Exact Solutions', 'FontSize', 16);
xlim([0, 1]);legend('Interpreter', 'latex', 'Location', 'best', 'FontSize', 12);set(gca, 'FontSize', 12, 'GridAlpha', 0.3, 'Box', 'on');
grid minor;

2.3、结果对比

可以看到数值解和解析解是几乎一致的：

3、全部代码

% Exercise 2.2 THE FINITE ELEMENT METHOD IN ELECTROMAGNETICS JIN JIAN MINGclc
clear
%% Solve
syms x c1 c2
% define the trail function and variational functional
phi = x+c1*sin(pi*x) +c2*sin(2*pi*x);
F=0.5*int((diff(phi, x))^2-2*(phi^2),x,0,1);
disp('F=');
pretty(simplify(F))% Find the extreme value of the functional. calculate the partial derivative of c
diff_c1 = diff(F, c1);  % ∂F/∂c1
diff_c2 = diff(F, c2);  % ∂F/∂c2% Solve the system of equations for variables [c1, c2]
equations = [diff_c1 == 0, diff_c2 == 0];
solutions = solve(equations, [c1, c2]);
disp('Solution:');
% Extract all solutions for c1 and c2
c1_val = double(solutions.c1);
c2_val = double(solutions.c2);
% Display results (both symbolic and numeric forms)
disp(['c1 = ', num2str(c1_val)]);
disp(['c2 = ', num2str(c2_val)]);phi_curve=subs(phi, {c1, c2}, {c1_val, c2_val});
phi_curve_exact=sin(sqrt(2)*x)/sin(sqrt(2));%% Draw figure
figure(1)
fplot(phi_curve, [0, 1], ...'Color', [0, 0.4470, 0.7410], ... % MATLAB经典蓝'LineStyle', '--', ...             % 虚线'LineWidth', 2.5, ...'Marker', 'o', ...                 % 圆形标记'MarkerSize', 7, ...'MarkerFaceColor', 'auto', ...     % 自动填充'MarkerEdgeColor', 'k', ...        % 黑色边框'DisplayName', sprintf('Numerical: $\\phi(x) = x + %.2f\\sin(\\pi x) + %.2f\\sin(2\\pi x)$', c1_val, c2_val));hold on
fplot(phi_curve_exact, [0, 1], ...'Color', [0.8500, 0.3250, 0.0980], ... % MATLAB经典红'LineStyle', '-', ...                   % 实线'LineWidth', 2.5, ...'Marker', 's', ...                      % 正方形标记'MarkerSize', 8, ...'MarkerFaceColor', 'auto', ...           % 自动填充'DisplayName', 'Exact: $\frac{\sin(\sqrt{2} x)}{\sin(\sqrt{2})}$');xlabel('Position (x)', 'FontSize', 14);
ylabel('Function Value \phi(x)', 'FontSize', 14);
title('Comparison of Numerical and Exact Solutions', 'FontSize', 16);
xlim([0, 1]);legend('Interpreter', 'latex', 'Location', 'best', 'FontSize', 12);set(gca, 'FontSize', 12, 'GridAlpha', 0.3, 'Box', 'on');
grid minor;

4、变分公式的证明