FDM.zip_fdm_打靶法_有限差分求解偏微分方程

format long clear clc a=0;b=1;Y0=1;h=1/121; t=a:h:b; n=(b-a)/h+1; Y(1)=Y0; %%R-K S(1)=1.5; S(2)=2.5; A=zeros(6,2); A(1,1)=a;A(1,2)=Y(1); z(1)=S(1); for i=2:n K(1,1)=f1(t(i-1),Y(i-1),z(i-1)); K(2,1)=f2(t(i-1),Y(i-1),z(i-1)); K(1,2)=f1(t(i-1)+h/2,Y(i-1)+K(1,1)*h/2,z(i-1)+K(2,1)*h/2); K(2,2)=f2(t(i-1)+h/2,Y(i-1)+K(1,1)*h/2,z(i-1)+K(2,1)*h/2); K(1,3)=f1(t(i-1)+h/2,Y(i-1)+K(1,2)*h/2,z(i-1)+K(2,2)*h/2); K(2,3)=f2(t(i-1)+h/2,Y(i-1)+K(1,2)*h/2,z(i-1)+K(2,2)*h/2); K(1,4)=f1(t(i-1)+h,Y(i-1)+K(1,3)*h,z(i-1)+K(2,3)*h); K(2,4)=f2(t(i-1)+h,Y(i-1)+K(1,3)*h,z(i-1)+K(2,3)*h); Y(i)=Y(i-1)+h/6*(K(1,1)+2*K(1,2)+2*K(1,3)+K(1,4)); z(i)=z(i-1)+h/6*(K(2,1)+2*K(2,2)+2*K(2,3)+K(2,4)); A(i,1)=t(i);A(i,2)=Y(i); end B(1)= Y(n); beta=1; e=0.5e-6; k=1; error=B(1)-beta; while abs(error)>=e k=k+1; z(1)=S(k); for i=2:n K(1,1)=f1(t(i-1),Y(i-1),z(i-1)); K(2,1)=f2(t(i-1),Y(i-1),z(i-1)); K(1,2)=f1(t(i-1)+h/2,Y(i-1)+K(1,1)*h/2,z(i-1)+K(2,1)*h/2); K(2,2)=f2(t(i-1)+h/2,Y(i-1)+K(1,1)*h/2,z(i-1)+K(2,1)*h/2); K(1,3)=f1(t(i-1)+h/2,Y(i-1)+K(1,2)*h/2,z(i-1)+K(2,2)*h/2); K(2,3)=f2(t(i-1)+h/2,Y(i-1)+K(1,2)*h/2,z(i-1)+K(2,2)*h/2); K(1,4)=f1(t(i-1)+h,Y(i-1)+K(1,3)*h,z(i-1)+K(2,3)*h); K(2,4)=f2(t(i-1)+h,Y(i-1)+K(1,3)*h,z(i-1)+K(2,3)*h); Y(i)=Y(i-1)+h/6*(K(1,1)+2*K(1,2)+2*K(1,3)+K(1,4)); z(i)=z(i-1)+h/6*(K(2,1)+2*K(2,2)+2*K(2,3)+K(2,4)); A(i,1)=t(i);A(i,2)=Y(i); end B(k)= Y(n); % if k<3 % S(k)=S(2); % else S(k+1)=S(k)-(B(k)-beta)/(B(k)-B(k-1))*(S(k)-S(k-1)); % end error=B(k)-beta; end Result.Sk=S(1:k); Result.YSk=B; Result ti=t'; Yi=Y'; Answer=table(ti,Yi) plot(t,Y,'r') hold on %% syms x; p=0.05; %p函数 r=x^2; %r函数 q=1; %q函数 f=0; h=(b-a)/n; %步长 value_of_f=zeros(n-1,1); %这是f diag_0=zeros(n-1,1); %确定A的对角元 diag_1=zeros(n-2,1); %确定A的偏离对角的上对角元 diag_2=zeros(n-2,1); X=a:h:b; a=1; b=1; for j=2:n %获取对角元素 diag_0(j-1)=((subs(p,{x},{(X(j+1)+X(j))/2}))+(subs(p,{x},{(X(j-1)+X(j))/2})))/(h^2)+(subs(q,{x},{X(j)})); end for j=3:n %获取A的第三条对角 diag_2(j-2)=-((subs(p,{x},{(X(j-1)+X(j))/2})))/(h^2)-subs(r,{x},{X(j)})/(2*h); end for j=2:n-1 %获取第二条对角 diag_1(j-1)=-((subs(p,{x},{(X(j-1)+X(j))/2})))/(h^2)+subs(r,{x},{X(j)})/(2*h); end for j=2:n value_of_f(j-1)=subs(f,{x},{X(j)}); end value_of_f(1)=value_of_f(1)+a*(subs(p,{x},{(X(2)+X(1))/2}))/(h^2); value_of_f(n-1)=value_of_f(n-1)+b*(subs(p,{x},{(X(n)+X(n+1))/2}))/(h^2); A=diag(diag_0)+diag(diag_1,1)+diag(diag_2,-1); u=A\value_of_f; U=[[a];u;[b]]; fprintf('%11.5f',U) plot(X,U,'b') legend('R-K','Diff')