一、顯示Euler
函數文件:Euler.m
1 function f=Euler(h,Y) 2 f(1,1)=Y(1)+h*(0.01-(1+(Y(1)+1000)*(Y(1)+1))*(0.01+Y(1)+Y(2))); 3 f(2,1)=Y(2)+h*(0.01-(1+Y(2)^2)*(0.01+Y(1)+Y(2)));
腳本文件:
1 tic; 2 clear 3 clc 4 %% 5 %顯示Euler方法求剛性微分方程,要求用Richardson外推法估計近似誤差從而控制步長 6 y(1:2,1)=[0;0];%初值 7 e=1e-5;%誤差過小 8 tol=1e-3;%指定的誤差 9 N=10;%節點的步數 10 h=1/N;%初始步長 11 t(1)=0; 12 i=1; 13 while t(i)+h<=1 14 k=1; 15 %%自動變步長 16 while k==1 17 y(1:2,i+1)=Euler(h,y(1:2,i));%符合誤差的數值解 18 y1_half=Euler(h/2,y(1:2,i));%半步長的中點數值解 19 y1_one=Euler(h/2,y1_half);%半步長的右端點的數值解 20 Estimate_error=2*norm(y(1:2,i+1)-y1_one);%中間估計誤差 21 if Estimate_error<tol%指定誤差 22 k=0;%步長相差不大,或者說正好在指定的誤差范圍內,則確定選擇h作為步長。 23 elseif Estimate_error<e%誤差過小 24 h=2*h; 25 else 26 h=h/2; 27 end 28 end 29 t(i+1)=t(i)+h; 30 i=i+1; 31 end 32 %% 33 %繪圖 34 plot(t,y,''); 35 xlabel('t'),ylabel('y(t) and z(t)'); 36 legend('y(t)','z(t)'); 37 title('Implicit Euler method for numerical solution of image'); 38 grid on; 39 toc;
效果圖:
二、隱式Euler:Euler.m
1 function X=Euler(t_h,u) 2 %隱式Euler(Newton迭代法) 3 %% 4 Tol=1e-5; 5 U=u; 6 x1=U-Jacobian(U,t_h)\F(U,u,t_h); 7 while (norm(x1-U,2)>=Tol) 8 %數值解的2范數是否在誤差范圍內 9 U=x1; 10 x1=U-Jacobian(U,t_h)\F(U,u,t_h); 11 end 12 X=x1;%不動點 13 %雅可比矩陣 14 function f=Jacobian(U,t_h) 15 f(1,1)=-t_h*((2*U(1)+1001)*(0.01+U(1)+U(2))+1+(U(1)+1000)*(U(1)+1))-1; 16 f(1,2)=-t_h*(1+(U(1)+1000)*(U(1)+1)); 17 f(2,1)=-t_h*(1+U(2)^2); 18 f(2,2)=-t_h*(2*U(2)*(0.01+U(1)+U(2))+(1+U(2)^2))-1; 19 20 %方程組 21 %% 22 function fun=F(U,u,t_h) 23 fun(1,1)=u(1)+t_h*(0.01-(1+(U(1)+1000)*(U(1)+1))*(0.01+U(1)+U(2)))-U(1); 24 fun(2,1)=u(2)+t_h*(0.01-(1+U(2)^2)*(0.01+U(1)+U(2)))-U(2);
腳本文件:
1 tic; 2 clear 3 clc 4 %隱式Euler方法求剛性微分方程,要求用Richardson外推法估計近似誤差從而控制步長 5 %% 6 y(1:2,1)=[0;0];%初值 7 e=1e-5;%誤差過小 8 tol=1e-3;%指定的誤差 9 N=100;%節點的步數 10 h=1/N;%初始步長 11 t(1)=0;%初始節點 12 i=1; 13 while t(i)+h<=1 14 k=1; 15 %自動變步長 16 while k==1 17 y(1:2,i+1)=Euler(h,y(1:2,i));%符合誤差的數值解 18 % y1_half=Euler(h/2,y(1:2,i));%半步長的中點數值解 19 y1_half=Euler(h/2,y(1:2,i));%半步長的右端點的數值解 20 Estimate_error=2*norm(y(1:2,i+1)-y1_half);%中間估計誤差 21 if Estimate_error<tol%指定誤差 22 k=0;%步長相差不大,或者說正好在指定的誤差范圍內,則確定選擇h作為步長。 23 elseif Estimate_error<e%誤差過小 24 h=2*h; 25 else%近似估計誤差大於指定誤差 26 h=h/2; 27 end 28 end 29 t(i+1)=t(i)+h; 30 i=i+1; 31 end 32 %繪圖 33 %% 34 plot(t,y); 35 xlabel('t'),ylabel('y(t) and z(t)'); 36 legend('y(t)','z(t)'); 37 title('Explicit Euler method for numerical solution of image'); 38 grid on ; 39 toc;
效果圖: