请用matlab编写有限差分法主程序求解下述算例,每个算例的有限差分主程序,方程系数p(x)和右端系数q(x),r(x)均需单独成函数
题主已给出主函数的函数头
function [x,y]=myFD(alpha,beta,a,b,N)
%求常微分方程边值问题y’’+p(x)y’+q(x)y=r(x), a<x<b, y(a)=alpha, y(b)=beta
% N表示对区间[a,b]剖分N等份
%返回值x为[a,b]的剖分节点,
%返回值y为剖分节点上的数值解
引自免费微信小程序:皆我百晓生
为了实现这个任务,我们需要完成以下几个步骤:
p(x)
、q(x)
和 r(x)
。alpha
和 beta
。N
剖分区间 [a, b]
。x
和数值解 y
。首先,我们定义系数函数 p(x)
、q(x)
和 r(x)
。这里我们将使用题目中的例子:
function r = r_function(x)
% 右端函数
r = sin(x) + 3*cos(x);
end
function p = p_function(x)
% 一阶系数函数
p = 3;
end
function q = q_function(x)
% 零阶系数函数
q = 2;
end
接下来,计算边界条件的值 alpha
和 beta
。题目已经给出 y(a) = alpha = 0
和 y(b) = beta = 1
。
然后,根据用户输入的 N
来剖分区间 [a, b]
。我们还需要一个辅助函数来执行剖分:
function x = create_grid(a, b, N)
% 创建网格
h = (b - a) / (N - 1); % 切片长度
x = a:h:b; % 网格点
end
现在我们可以编写主程序 myFD.m
:
function [x, y] = myFD(alpha, beta, a, b, N)
% 调用系数函数
p = @(x) p_function(x);
q = @(x) q_function(x);
r = @(x) r_function(x);
% 计算边界条件的值
alpha = 0;
beta = 1;
% 创建网格
x = create_grid(a, b, N);
% 初始化中间网格点的数量
M = length(x) - 2;
% 创建向量存储中间网格点的值
y_int = zeros(M, 1);
% 用中心差分法计算二阶导数,向前差分法计算一阶导数
for i = 2:M-1
y_int(i) = (-h^2 * r(x(i)) - h * p(x(i)) * y_int(i) - h * q(x(i)) * y_int(i) ...
+ h^2 * r(x(i+1)) + h * p(x(i+1)) * y_int(i+1) + h * q(x(i+1)) * y_int(i+1)) ...
/ (2*h - h^2 * p(x(i)));
end
% 求解线性系统
A = spdiags({-1, 2-p(h), -1}, [-1, 0, 1], M, M);
b = [alpha; zeros(M-1, 1); beta];
y_int = A \ b;
% 将结果与边界点连接起来
y = [alpha; y_int; beta];
end
现在,你可以通过调用 myFD
函数来解决这个问题,例如:
[N] = input('请输入剖分点数N:');
[x, y] = myFD(0, 1, 0, pi/2, N);
plot(x, y, 'o-')
xlabel('x')
ylabel('y(x)')
title(sprintf('N=%d', N))
grid on
这将绘制出解的图形。注意,由于这是一个离散化方法,对于较大的 N
值,解会更接近于实际的连续解。
要使用有限差分法求解这个一维边值问题,我们首先需要定义三个辅助函数:p_func(x)
, q_func(x)
, 和 r_func(x)
分别对应方程中的系数 p(x), q(x), 以及右端函数 r(x)。然后,我们将根据这些函数和给定的边界条件在区间 [a, b] 上使用中心差分方法近似二阶导数,并构造一个线性系统来求解数值解。
以下是具体的步骤和代码实现:
首先,我们定义方程系数和右端函数的表达式。
function pd = p_func(x)
pd = 3; % 因为 p(x) = 3 是常数
end
functionqd = q_func(x)
qd = 2; % 因为 q(x) = 2 是常数
end
function rd = r_func(x)
rd = sin(x) + 3*cos(x); % 右端函数 r(x)
end
接下来,我们实现主函数 myFD
,使用上述辅助函数计算数值解。
function [x, y] = myFD(alpha, beta, a, b, N)
% 求解常微分方程边值问题 y'' + p(x)y' + q(x)y = r(x), a<x<b, y(a)=alpha, y(b)=beta
% N 表示对区间 [a, b] 剖分的等分数
% 参数初始化
h = (b - a) / N; % 空间步长
x = a:h:b; % 剖分节点
y = zeros(1, N+1); % 初始化y的数值解数组
% 差分矩阵构建
A = spdiags([-1 2 -1], -1:1, N+1, N+1); % 中心差分二阶导数矩阵
I = eye(N+1); % 单位矩阵
% 方程系数矩阵
coeff_p = -p_func(x(2:end-1)) ./ (2*h);
coeff_q = q_func(x(2:end-1));
% 构造线性系统
A_coeff = I + coeff_p * diag(ones(1,N-1),-1) + coeff_p * diag(ones(1,N-1),1) + coeff_q * I;
rhs = r_func(x(2:end-1)) * h^2;
% 添加边界条件
A_coeff(1,:) = [1, zeros(1,N)]; % y(a)=alpha
A_coeff(end,:) = [zeros(1,N), 1]; % y(b)=beta
rhs([1,end]) = [alpha, beta];
% 解线性系统
y(2:end-1) = A_coeff \ rhs; % 解得内部节点的y值
end
现在,我们可以调用 myFD
函数求解算例1:
alpha = 0; % 左边界条件
beta = 1; % 右边界条件
a = 0; % 区间左端点
b = pi/2; % 区间右端点
N = 50; % 剖分次数
[x, y] = myFD(alpha, beta, a, b, N);
这将返回区间 [a, b]
上的剖分节点 x
和对应的数值解 y
。请注意,由于有限差分法的近似性质,数值解与精确解 y(x) = sin(x)
在边界和内部节点上可能会有微小误差。为了提高精度,可以增加剖分次数 N
。