用差分法求出数值解之后,求与精确解误差出问题了。

2 次查看(过去 30 天)
求误差得到的结果是原函数的图不知道那里的语句出问题了,图片左为数值解图像右边误差除了问题,求解
程序如下:
clear all; close all;
a=4; h=0.05; x=[0:h:1];
tau=0.0125; t=[0:tau:1];
s=a*tau/h;
N=length(x)-1; M=length(t)-1;
[T X]=meshgrid(t,x);
%构造矩阵
e=s^2*ones(N-1,1);
A=spdiags([e 2*(1-e) e],[-1 0 1],N-1,N-1);
%设置初始条件和边界条件
u=zeros(N+1,M+1);
u(:,1)=sin(pi*x); u(:,2)=0;
u(1,:)=0; u(end,:)=0;
%有限差分
for n=2:M
u(2:N,n+1)=A*u(2:N,n)-u(2:N,n-1);
u(2,n+1)=u(2,n+1)+s^2*u(1,n);
u(N,n+1)=u(N,n+1)+s^2*u(end,n);
end
%画图
subplot(2,2,1)
mesh(t,x,u), view(10,10), xlabel t, ylabel x, zlabel u
subplot(2,2,2)
%与解析解的误差
m=sin(pi*X).*cos(4*pi*T)
mesh(t,x,u-m), view(10,8), axis([0 1 0 1 -10 10])
xlabel t, ylabel x, zlabel Error

采纳的回答

华纳娱乐游戏网址【359663.tv】
直接写了个对的,也没写成复杂的矩阵运算格式,你去看看我推荐的讲偏微分方程的书籍的对应章节再对照代码吧。
clear; clc; close all;
c = 4;
X_min = 0; X_Step = 1/10; X_max = 1;
T_min = 0; T_Step = 0.025; T_max = 1;
F = @( x ) sin( pi * x );
G = @( x ) 0 * x;
X_Vector = X_min : X_Step : X_max;
T_Vector = T_min : T_Step : T_max;
lambda = c * T_Step / X_Step;
U_Matrix = NaN * ones( numel( X_Vector ), numel( T_Vector ) );
U_Matrix( 1, : ) = 0; U_Matrix( end, : ) = 0;
U_Matrix( :, 1 ) = F( X_Vector ).';
U_Matrix( 2 : end - 1, 2 ) = ( lambda^2 * ( F( X_Vector( 1 : end - 2 ) ) + F( X_Vector( 3 : end ) ) ) / 2 + ...
    ( 1 - lambda^2 ) * F( X_Vector( 2 : end - 1 ) )  + T_Step * G( X_Vector( 2 : end - 1 ) ) ).';
for jj = 3 : 1 : numel( T_Vector )
    U_Matrix( 2 : end - 1, jj ) = lambda^2 * ( U_Matrix( 1 : end - 2, jj - 1 ) + U_Matrix( 3 : end, jj - 1 ) ) + ...
        ( 1 - lambda^2 ) * 2 * U_Matrix( 2 : end - 1, jj - 1 )  - U_Matrix( 2 : end - 1, jj - 2 );
end
[ X_Matrix, T_Matrix ] = meshgrid( X_Vector, T_Vector );
U_Analytic = @( x, t ) sin( pi * x ) .* cos( c * pi * t );
U_AnaMatrix = U_Analytic( X_Matrix, T_Matrix );
%%
figure;
subplot( 1, 2, 1 )
surf( X_Matrix, T_Matrix, U_Matrix.',...
    'FaceColor', 'r', 'FaceAlpha', 0.7, 'MeshStyle', 'default', 'EdgeColor', 'k', 'EdgeAlpha', 0.6, ...
    'AlignVertexCenters', 'on', 'LineStyle', ':', 'LineWidth', 0.8 );
xlabel( 'x', 'FontSize', 18 ); ylabel( 't', 'FontSize', 18 ); zlabel( 'u_{numeric}', 'FontSize', 18 );
axis equal; grid on;
subplot( 1, 2 ,2 )
surf( X_Matrix, T_Matrix, U_AnaMatrix,...
    'FaceColor', 'b', 'FaceAlpha', 0.7, 'MeshStyle', 'default', 'EdgeColor', 'k', 'EdgeAlpha', 0.6, ...
    'AlignVertexCenters', 'on', 'LineStyle', ':', 'LineWidth', 0.8 );
xlabel( 'x', 'FontSize', 18 ); ylabel( 't', 'FontSize', 18 ); zlabel( 'u_{analytic}', 'FontSize', 18 );
axis equal; grid on;

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Logical 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!