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

求误差得到的结果是原函数的图不知道那里的语句出问题了,图片左为数值解图像右边误差除了问题,求解
程序如下:
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

 採用された回答

0 投票

直接写了个对的,也没写成复杂的矩阵运算格式,你去看看我推荐的讲偏微分方程的书籍的对应章节再对照代码吧。
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 件)

カテゴリ

ヘルプ センター および File ExchangePhased Array System Toolbox についてさらに検索

タグ

Community Treasure Hunt

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

Start Hunting!