Matlab:顯(隱)式Euler和Richardson外推法變步長求解剛性問題


一、顯示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;

效果圖:

 


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM