crank_nicolson3_12mended2.m

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

function [e]=Crank_Nicolson3_12mended2(h,tao)
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
k=1/tao+1;
  for i=1:1/h+1
      W(k,i)=exp(X(1,i))*sin(0.5-T(1,k));E(1,i)=abs(W(k,i)-U(k,i));
  end
  plot(X,E);
  hold on;