homework of numerical solution of PDEs

来源:互联网 发布:分析数据的统计处理 编辑:程序博客网 时间:2024/05/22 13:23

Problem set

find the numerical solution of the following partial differential equation:


1.C - N  Scheme


%二阶波动方程的CN格式A=[0,-1;-1,0;];I=[1,0;0,1];N=100;%空间划分多少步h=1.0/100;t=0.0001;r=t/h;step=200;%时间步数U=zeros(2,N+1);d=zeros(2,N+1);a11=1:N+1;a12=1:N+1;a21=1:N+1;a22=1:N+1;b11=1:N+1;b12=1:N+1;b21=1:N+1;b22=1:N+1;c11=1:N+1;c12=1:N+1;c21=1:N+1;c22=1:N+1;for i=1:1:N+1    a11(i)=-0.25*A(1,1)*r;b11(i)=I(1,1);c11(i)=0.25*A(1,1)*r;    a12(i)=-0.25*A(1,2)*r;b12(i)=I(1,2);c12(i)=0.25*A(1,2)*r;    a21(i)=-0.25*A(2,1)*r;b21(i)=I(2,1);c21(i)=0.25*A(2,1)*r;    a22(i)=-0.25*A(2,2)*r;b22(i)=I(2,2);c22(i)=0.25*A(2,2)*r;endV=zeros(N+1,step+1);W=zeros(N+1,step+1);u=zeros(N+1,step+1);%实际最终求的解%第一步,求解得到V,M%初始化,进行赋值x=0:h:1;V(:,1)=-exp(x);%初始值W(:,1)=exp(x);%初始值for k=1:1:step    %边界值    V(1,k)=-exp(-(k-1)*t);    V(N+1,k)=-exp(1-(k-1)*t);    %将V赋给U(1,).将W赋给U(2,)    U(1,:)=V(:,k);    U(2,:)=W(:,k);    %追赶法求解        for j=2:1:N        d(:,j)=U(:,j)+0.25*A*r*(U(:,j+1)-U(:,j-1));    end    %边界条件    d(:,1)=U(:,1)+0.25*A*r*U(:,2);    d(:,N+1)=U(:,N+1)-0.25*A*r*U(:,N);    for j=2:1:N+1        A1=[a11(j),a12(j);a21(j),a22(j);];        B=[b11(j-1),b12(j-1);b21(j-1),b22(j-1);];        A1=B\A1;        a11(j)=A1(1,1);a12(j)=A1(1,2);a21(j)=A1(2,1);a22(j)=A1(2,2);        C=[c11(j-1),c12(j-1);c21(j-1),c22(j-1);];        B=[b11(j),b12(j);b21(j),b22(j);];        B=B-C*A1;        b11(j)=B(1,1);b12(j)=B(1,2);b21(j)=B(2,1);b22(j)=B(2,2);        d(:,j)=d(:,j)-A1*d(:,j-1);    end    B=[b11(N+1),b12(N+1);b21(N+1),b22(N+1);];    d(:,N+1)=B\d(:,N+1);    for j=N:-1:1        C=[c11(j),c12(j);c21(j),c22(j);];        B=[b11(j),b12(j);b21(j),b22(j);];        d(:,j)=B\(d(:,j)-C*d(:,j+1));    end     V(:,k+1)=d(1,:);%求得VW    W(:,k+1)=d(2,:);end%第二步从VW求得uu(:,1)=exp(x);for k=1:1:step    u(1,k)=exp(-(k-1)*t);    u(N+1,k)=exp(1-(k-1)*t);    for i=1:1:N        u(i+1,k+1)=u(i,k)+t*V(i,k)+h*W(i,k);     end    u(1,k+1)=exp(-(k)*t);    u(N+1,k+1)=exp(1-(k)*t);end
2.C-F-L Scheme


%CFL格式求解N=100;%节点数M=200;%时间步V=zeros(N+1,M+1);W=zeros(N+1,M+1);U=zeros(N+1,M+1);t=0.0001;h=1.0/N;r=t/h;v=zeros(N+1,1);w=zeros(N+1,1);%初始化数据x=0:h:1.0;V(:,1)=-exp(x);x=0.5*h:h:1+h;W(:,1)=exp(x);%求解V Wfor k=1:1:M    V(1,k)=-exp(-(k-1)*t);    V(N+1,k)=-exp(1-(k-1)*t);    v=V(:,k);    w=W(:,k);    %求v n+1    for i=2:1:N+1        V(i,k+1)=v(i)+r*(w(i)-w(i-1));    end    V(1,k+1)=-exp(-k*t);    V(1,k+1)=-exp(1-k*t);    %求W n+1    for i=2:1:N+1        W(i,k+1)=w(i)+r*(V(i,k+1)-V(i-1,k+1));    end        W(1,k+1)=w(1)+r/2*(t+t);%边界值endU(:,1)=exp(x);for k=1:1:M    U(1,k)=exp(-(k-1)*t);    U(N+1,k)=exp(1-(k-1)*t);    for i=1:1:N        U(i+1,k+1)=U(i,k)+t*V(i,k)+h*W(i,k);    end    U(1,k+1)=exp(-(k)*t);    U(N+1,k+1)=exp(1-(k)*t);end
3 隐格式
%Yin Scheme Second Order EquN=100;%节点数M=200;%时间步h=1.0/N;t=0.0001;A=1.0;r=A*t/h;V=zeros(N+1,M+1);W=zeros(N+1,M+1);v=zeros(N+1,1);w=zeros(N+1,1);u=zeros(N+1,M+1);a=zeros(N+1,1);b=zeros(N+1,1);c=zeros(N+1,1);d=zeros(N+1,1);for i=1:1:N+1    a(i)=-A*A*r*r/4;    b(i)=1+A*A*r*r/2;    c(i)=-A*A*r*r/4;end%初始值x=0:h:1;V(:,1)=-exp(x);m=0.5*h;x=m:h:1+h;W(:,1)=exp(x);%求解for k=1:1:M    V(1,k)=-exp(-(k-1)*t);    V(N+1,k)=-exp(1-(k-1)*t);    v=V(:,k);    w=W(:,k);    %先求v,n+1    for i=2:1:N        d(i)=v(i)+A*r*(w(i+1)-w(i))+A*A*r*r/4*(v(i+1)-2*v(i)+v(i-1));    end    d(1)=v(1)+A*r*(w(2)-w(1))+A*A*r*r/4*(v(2)-2*v(1));    d(N+1)=v(N+1)+A*r*(-w(N+1))+A*A*r*r/4*(-2*v(N+1)+v(N));    for i=2:N+1        a(i)=a(i)/b(i-1);        b(i)=b(i)-c(i-1)*a(i);        d(i)=d(i)-a(i)*d(i-1);    end    d(N+1)=d(N+1)/b(N+1);    for i=N:-1:1        d(i)=(d(i)-c(i)*d(i+1))/b(i);    end    V(:,k+1)=d;    V(1,k+1)=-exp(-(k)*t);    V(N+1,k+1)=-exp(1-(k)*t);    d=V(:,k+1);    %求W n+1    for i=2:1:N+1        W(i,k+1)=w(i)+A*r/2*(d(i)-d(i-1)+v(i)-v(i-1));    end    W(1,k+1)=w(1)+A*r/2*(t+t);    end%求uu(:,1)=exp(x);for k=1:1:M    u(1,k)=exp(-(k-1)*t);    u(N+1,k)=exp(1-(k-1)*t);    for i=1:1:N        u(i+1,k+1)=u(i,k)+t*V(i,k)+h*W(i,k);    end    u(1,k+1)=exp(-(k)*t);    u(N+1,k+1)=exp(1-(k)*t);end




0 0