matlab-programs.rar_numerical analysis

clear all close all format long L=1;N=80;h=L/N;tt=5/10;l=1/120; A=eye(N,N); j=1; A(1,j)=-25;A(1,j+1)=48;A(1,j+2)=-36;A(1,j+3)=16;A(1,j+4)=-3; A(2,j)=-18;A(2,j+1)=20;A(2,j+2)=-12;A(2,j+3)=9;A(2,j+4)=-2; A(3,j)=-9;A(3,j+1)=12;A(3,j+2)=-20;A(3,j+3)=18;A(3,j+4)=-3; A(4,j)=-16;A(4,j+1)=36;A(4,j+2)=-48;A(4,j+3)=25; for r=5:N; A(r,r-4)=3;A(r,r-3)=-16;A(r,r-2)=36;A(r,r-1)=-48;A(r,r)=25; end A=(1/(12*h))*A; eig(A); a1=64/25;a2=7/3;a3=547/600; r1=-.937580908085238; r2=-1.81471985800593; r3=-2.00000000000000; r4=-2.26051974672921; col_vec1=[]; for x=h:h:L; col_vec1=[col_vec1;exp(x)]; end gx=col_vec1; I=eye(N,N); X1=I-(l/r1)*A; X2=I-(l/r2)*A; X3=I-(l/r3)*A; X4=I-(l/r4)*A; num1=1-(1-a1)*r1+(.5-a1+a2)*(r1)^2-(1/6-.5*a1+a2-a3)*(r1)^3; num2=1-(1-a1)*r2+(.5-a1+a2)*(r2)^2-(1/6-.5*a1+a2-a3)*(r2)^3; num3=1-(1-a1)*r3+(.5-a1+a2)*(r3)^2-(1/6-.5*a1+a2-a3)*(r3)^3; num4=1-(1-a1)*r4+(.5-a1+a2)*(r4)^2-(1/6-.5*a1+a2-a3)*(r4)^3; den1=(1-r1/r2)*(1-r1/r3)*(1-r1/r4); den2=(1-r2/r1)*(1-r2/r3)*(1-r2/r4); den3=(1-r3/r1)*(1-r3/r2)*(1-r3/r4); den4=(1-r4/r1)*(1-r4/r2)*(1-r4/r3); p1=num1/den1; p2=num2/den2; p3=num3/den3; p4=num4/den4; num5=3-(-19+78*a1-216*a2+324*a3)*r1+(3-8*a1+12*a2)*(r1)^2; num6=3-(-19+78*a1-216*a2+324*a3)*r2+(3-8*a1+12*a2)*(r2)^2; num7=3-(-19+78*a1-216*a2+324*a3)*r3+(3-8*a1+12*a2)*(r3)^2; num8=3-(-19+78*a1-216*a2+324*a3)*r4+(3-8*a1+12*a2)*(r4)^2; p5=num5/den1; p6=num6/den2; p7=num7/den3; p8=num8/den4; num9=2-(16-56*a1+144*a2-216*a3)*r1+(1-4*a1+12*a2-24*a3)*(r1)^2; num10=2-(16-56*a1+144*a2-216*a3)*r2+(1-4*a1+12*a2-24*a3)*(r2)^2; num11=2-(16-56*a1+144*a2-216*a3)*r3+(1-4*a1+12*a2-24*a3)*(r3)^2; num12=2-(16-56*a1+144*a2-216*a3)*r4+(1-4*a1+12*a2-24*a3)*(r4)^2; p9=num9/den1; p10=num10/den2; p11=num11/den3; p12=num12/den4; num13=1-(-7+26*a1-72*a2+108*a3)*r1+(-1+4*a1-12*a2+24*a3)*(r1)^2; num14=1-(-7+26*a1-72*a2+108*a3)*r2+(-1+4*a1-12*a2+24*a3)*(r2)^2; num15=1-(-7+26*a1-72*a2+108*a3)*r3+(-1+4*a1-12*a2+24*a3)*(r3)^2; num16=1-(-7+26*a1-72*a2+108*a3)*r4+(-1+4*a1-12*a2+24*a3)*(r4)^2; p13=num13/den1; p14=num14/den2; p15=num15/den3; p16=num16/den4; num17=6-(44-168*a1+432*a2-648*a3)*r1+(11-44*a1+132*a2-216*a3)*(r1)^2-(2-8*a1+24*a2-48*a3)*(r1)^3; num18=6-(44-168*a1+432*a2-648*a3)*r2+(11-44*a1+132*a2-216*a3)*(r2)^2-(2-8*a1+24*a2-48*a3)*(r2)^3; num19=6-(44-168*a1+432*a2-648*a3)*r3+(11-44*a1+132*a2-216*a3)*(r3)^2-(2-8*a1+24*a2-48*a3)*(r3)^3; num20=6-(44-168*a1+432*a2-648*a3)*r4+(11-44*a1+132*a2-216*a3)*(r4)^2-(2-8*a1+24*a2-48*a3)*(r4)^3; p17=num17/den1; p18=num18/den2; p19=num19/den3; p20=num20/den4; col_vec2=[]; col_vec2=[col_vec2;0]; col_vec2=[col_vec2;-1/(4*h)]; col_vec2=[col_vec2;-1/(6*h)]; col_vec2=[col_vec2;-1/(4*h)]; for j=5:N; col_vec2=[col_vec2;0]; j=j+1; end for t=0:1:tt-l; t1=t; v1=col_vec2*exp(-t1); t2=t1+l/3; v2=col_vec2*exp(-t2); t3=t2+l/3; v3=col_vec2*exp(-t3); t4=t3+l/3; v4=col_vec2*exp(-t4); w1=(l/48)*(2*p5*v1+9*p9*v2+18*p13*v3+p17*v4); w2=(l/48)*(2*p6*v1+9*p10*v2+18*p14*v3+p18*v4); w3=(l/48)*(2*p7*v1+9*p11*v2+18*p15*v3+p19*v4); w4=(l/48)*(2*p8*v1+9*p12*v2+18*p16*v3+p20*v4); if t==0; ut=gx; else ut=un; end y1=inv(X1)*(p1*ut+w1); y2=inv(X2)*(p2*ut+w2); y3=inv(X3)*(p3*ut+w3); y4=inv(X4)*(p4*ut+w4); un=y1+y2+y3+y4; end unum=un uext=gx*exp(-tt); err=uext-unum; col_vecnum=[]; col_vecnum=[col_vecnum;unum]; y=col_vecnum; col_vecext=[]; col_vecext=[col_vecext;uext]; x=linspace(0,1,N); w=col_vecext; plot(x,y,'-');xlabel('x');ylabel('u'); plot(x,w,'--');xlabel('x');ylabel('u'); plot(x,y,'-',x,w,'--');