crank_nicolson3_12.m

利用MATLAB解决偏微分方程中的一些问题,如插值求解

function[Y1,Y2,E1]=Crank_Nicolson3_12(h,tao,A)
a=2;r=a*tao/h^2;p=r/2;
for k=1:1/tao+1
  for i=1:1/h+1
      U(k,i)=0;
  end
end
for i=1:1/h+1
    X(1,i)=(i-1)*h;U(1,i)=exp(X(1,i))*sin(0.5);
end

for k=1:1/tao+1
    T(1,k)=(k-1)*tao;U(k,1)=sin(0.5-T(1,k));U(k,1/h+1)=exp(1)*sin(0.5-T(1,k));
end
for i=1:1/h-1
    for j=1:1/h-1
        Z(i,j)=0; Z(i,j)=0;
    end
end
for i=1:1/h-1
    Z(i,i)=1+2*p;W(i,i)=1-2*p;
end
for i=2:1/h-1
    Z(i,i-1)=-p;W(i,i-1)=p;
    Z(i-1,i)=-p; W(i-1,i)=p;
end
    for k=2:1/tao+1
    for j=2:1/h
        B(j-1,1)=tao*(-exp(X(1,j))*(cos(0.5-T(1,k-1))+2*sin(0.5-(T(1,k-1)+T(1,k))/2)));
end
    B(1,1)=B(1,1)+p*(U(k-1,1)+U(k,1));
    B(1/h-1,1)=B(1/h-1,1)+p*(U(k-1,1/h+1)+U(k,1/h+1));
    U([k],2:1/h)=(inv(Z)*((W*(U([k-1],2:1/h))')+B))';
end
for j=1:length(A)
Y1(1,j)=U(ceil((A(1,j)/tao)+1),ceil((0.5/h)+1));Y2(1,j)=exp(0.5)*sin(0.5-A(1,j));E1(1,j)=abs(Y1(1,j)-Y2(1,j));
end