优·matlab程序范例.docVIP

1.二阶椭圆方程 function eqsolv （1）.% -dd u/dxx + u=f on [aa,bb] with u(a)=sin(aa),u(b)=sin(bb) format long %format short aa=-1; bb=1; h=0.00005; N=(bb-aa)/h x=zeros(N+1,1); RHS=zeros(N-1,1); un=zeros(N+1,1); ue=zeros(N+1,1); for k=1:N+1 x(k)=aa+h*(k-1); ue(k)=sin(x(k)); end for k=1:N-1 RHS(k)=2*h*2*sin(x(k+1)); end RHS(1)=RHS(1)+ue(1)*2/h; RHS(N-1)=RHS(N-1)+ue(N+1)*2/h; % to use gauss method %A=zeros(N-1,N); %for k=1:N-2 % A(k,k+1)=-2/h; % A(k,k)=4/h+2*h; % A(k+1,k)=-2/h; %end % A(N-1,N-1)=4/h+2*h; %A(1:N-1,N)=RHS; %un(2:N)=gaussc(A); A=zeros(N-1,3); for k=1:N-1 A(k,2)=4/h+2*h; % A(k,k)=4/h+2*h; A(k,1)=-2/h; % A(k+1,k)=-2/h; A(k,3)=-2/h; % A(k,k+1)=-2/h; end L=zeros(N-1,1); R=zeros(N-1,2); R(1,1)=A(1,2); % R(1,1)=A(1,1); for k=1:N-2 R(k,2)=A(k,3); % R(k,2)=A(k,k+1); end for k=2:N-1 L(k)=A(k-1,1)/R(k-1,1); % L(k)=A(k,k-1)/R(k-1,1); R(k,1)=A(k,2)-L(k)*R(k-1,2); % R(k,1)=A(k,k)-L(k)*R(k-1,2); end % to solve LY=d y=zeros(N-1,1); y(1)=RHS(1); for k=2:N-1 y(k)=RHS(k)-L(k)*y(k-1); end % % to solve UX=Y un(N)=y(N-1)/R(N-1,1); for k=N-1:-1:2 un(k)=(y(k-1)-R(k-1,2)*un(k+1))/R(k-1,1); end % to compute the max err un(1)=ue(1); un(N+1)=ue(N+1); un=un-ue; err_max=abs(un(1)); for k=2:N+1 if err_maxabs(un(k)) err_max=abs(un(k)); end end err_max （2）.function ejtyfc % -d(p*du/dx)/dx+q*u=f on [aa,bb] with u(a)=sin(aa),u(b)=sin(bb) % p(a)=cos(aa),p(b)=cos(bb) q(a)=aa^2,q(b)=bb^2 format long format short aa=-1; bb=1; h=0.0005; N=(bb-aa)/h; x=zeros(N+1,1); x1=zeros(N,1); g=zeros(N-1,1); ue=zeros(N+1,1); un=zeros(N+1,1); for i=1:N+1 x(i)=aa+h*(i-1); ue(i)=sin(x(i)); end for i=1:N x1(i)=(x(i)+x(i+1))/2; end for i=1:N-1 g(i)=2*h*(sin(2*x(i+1))+x(i+1)^2*sin(x(i+1))); end g(1)=g(1)+ue(1)*2*cos(x1(1))/h; g(N-1)=g(N-1)+ue(N+1)*2*cos(x1(N))/h; if(0) A=zeros(N-1,N); for i=1:N-2 A(i,i)=2*(cos(x1(i+1))+cos(x1(i)))/h+2*h*x(