difds2.m

用LAPLACE方法求解多维沃尔泰拉方程

% matlab programm to compute partitial derivative equation with the difference method
% date: 2003-3-25,time:8.00pm
% equation:du/dt-integal[(t-s)^2*d^2u/dx^2,{0,t}]=sin(tx)
% bound: u(0,t)=u(1,t)=0
% initial: u(x,0)=x(x-1)
% 0<x<1,0<=t<=infinit
clear all;
%global J ux0 uxN;
t0=0;
t1=1;
x0=0;
x1=1;
J=input('Please input the x grid numer:\n')
N=input('Please input the t grid numer:\n')
dt=(t1-t0)/N;
dx=(x1-x0)/J;
lambda=dt/dx;
t=[t0:dt:t1];
x=[x0:dx:x1];
fx=[x0+dx:dx:x1-dx];
u=zeros(N+1,J+1);
u(1,:)=x.*(x-ones(1,J+1));
u(:,1)=zeros(N+1,1);
u(:,J+1)=zeros(N+1,1);
A=eye(J-1)*(1+2*lambda^2/sqrt(dt))+sparse(2:J-1,1:J-2,ones(1,J-2),J-1,J-1)*(-lambda^2/sqrt(dt))+(sparse(2:J-1,1:J-2,ones(1,J-2),J-1,J-1)*(-lambda^2/sqrt(dt)))';
IA=inv(A);
B=eye(J-1)*(-2)+sparse(2:J-1,1:J-2,ones(1,J-2),J-1,J-1)+(sparse(2:J-1,1:J-2,ones(1,J-2),J-1,J-1))';
tic;
for i=2:N+1
  tt=sin(fx'*(i-1)*dt)*dt;
  tt(1)=tt(1)+u(i-1,1)*lambda^2/sqrt(dt);
  tt(J-1)=tt(J-1)+u(i-1,J+1)*lambda^2/sqrt(dt);
  u0=u(i-1,2:J)';
  ut=zeros(J-1,1);
  for j=2:i-1
      u1=u(j,2:J)';
      uf=(B*u1)*lambda^2/sqrt((i-j+1)*dt);
      ut=ut+uf;
  end    
  u0=u0+tt;
  ut=ut+u0;
  u1=IA*ut;
  u(i,2:J)=u1';
end
u
%up=(vpa(u,20))
t0=toc;
t1=fix(clock);
clear ut0 ux0 uxN A B u0 u1 uf ut tt
figure(1);
title(['\The solution surface']);
mesh(u);grid on;
fid=fopen('difds2.txt','w+t');
fprintf(fid,'%8f\n',t0);
fprintf(fid,'%20d\n',t1);
fprintf(fid,'NxJ=%-5dx%5d\n',N,J);
fprintf(fid,'{');
for i=1:N+1
    fprintf(fid,'{');
    for j=1:J+1    
       fprintf(fid,'%20.8e,',u(i,j));
    end   
    fprintf(fid,'},\n');
end
fprintf(fid,'}');
%fprintf(fid,'%-20.8e\n',u);
fclose(fid);