%%??p281頁?對角隱式辛積分R-K方法?精度為4階??1e-4
%??調用格式：[TZ]?=?DiagImpSym（functspanh0z0ers）
%??其中：????func為微分方程，i.e.?func?=?@（tz）cos（t）
%????????????tspan為起止時間，tspan=[0?2*pi]?h0為步長
%????????????z0為微分方程狀態變量的初始值z0?=?0?ers為迭代精度誤差上限
function?[TZ]?=?DiagImpSym（functspanh0z0ers）
???%%?公式系數
%????b2?=?-0.53652708;
%????b3?=?2.37893931;
%????b4?=?1.8606818856;
%????b1?=?1.0-b2-b3-b4;
???load?DiagImpSymCoef.mat
???b1?=?X（1）;
???b2?=?X（2）;
???b3?=?X（3）;
???b4?=?X（4）;
???z0?=?reshape（z0length（z0）1）;
???%%?時間節點
???t0?=?tspan（1）;
???tF?=?tspan（end）;
???if（t0>=tF）
???????error（‘起止時間輸入錯誤！‘）;
???end
%????h0?=?1e-5;%步長
???if（h0>tF-t0）
???????error（‘輸入的步長過大！‘）
???end
???nh?=?floor（（tF-t0）/h0）;
???if（tF>=t0+nh*h0）
??????T?=?linspace（t0tFnh+1）.‘;
??????h?=?（tF-t0）/（nh+1）;
???else
??????T?=?linspace（t0tFnh）.‘;
??????h?=?（tF-t0）/nh;
???end
???num_T?=?length（T）;
??
???%%?微分方程的輸入個數
???num_var?=?nargin（func）;
???if（num_var~=2）???????
??????error（‘微分方程的輸入個數錯誤！‘）;?
???end
???%%?求解微分方程
???zk?=?z0;
???Z?=?zeros（num_Tlength（z0））;
???Z（1:）?=?zk;
???for?k?=?1:num_T-1???
???????%?t
???????tk?=?T（k）;
???????tk1?=?T（k+1）;
???????h?=?tk1-tk;
???????%?Y1
???????Y10?=?zk+b1*h/2*feval（functk+b1*h/2zk）;%時間節點的問題，初值問題
???????while（1）
???????????Y1?=?zk+b1*h/2*feval（functk+b1*h/2Y10）;
???????????if（norm（Y1-Y10）/norm（Y1）???????????????break;
???????????else
???????????????Y10?=?Y1;?
???????????end???????????
???????end
???????%?Y2
????????Y20?=?2*Y1-zk+b2*h/2*feval（functk+（2*b1+b2）*h/22*Y1-zk）;%時間節點的問題，初值問題
???????while（1）
???????????Y2?=?2*Y1-zk+b2*h/2*feval（functk+（2*b1+b2）*h/2Y20）;
???????????if（norm（Y2-Y20）/norm（Y2）???????????????break;
???????????else
???????????????Y20?=?Y2;?
???????????end???????????
???????end
???????%?Y3
????????Y30?=?2*Y2-（2*Y1-zk）+b3*h/2*feval（functk+（2*b1+2*b2+b3）*h/22*Y2-（2*Y1-zk））;%時間節點的問題，初值問題
???????while（1）
????????????Y3?=?2*Y2-（2*Y1-zk）+b3*h/2*feval（functk+（2*b1+2*b2+b3）*h/2Y30）;
???????????if（norm（Y3-Y30）/norm（Y3）???????????????break;
???????????else
???????????????Y30?=?Y3;?
???????????end???????????
???????end
???????%?Y4
????????Y40?=?2*Y3-（2*Y2-2*Y1+zk）+b4*h/2*feval（functk+（2*b1+2*b2+2*b3+b4）*h/22*Y3-（2*Y2-2*Y1+zk））;%時間節點的問題，初值問題
???????while（1）
????????????Y4?=?2*Y3-（2*Y2-2*Y1+zk）+b4*h/2*feval（functk+（2*b1+2*b2+2*b3+b4）*h/2Y40）;
???????????if（norm（Y4-Y40）/norm（Y4）???????????????break;
???????????else
???????????????Y40?=?Y4;?
???????????end???????????
???????end
???????%?zk1
???????zk1?=?2*Y4-（2*Y3-2*Y2+2*Y1-zk）;
???????Z（k+1:）?=?zk1;
???????zk?=?zk1;
???end
end

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----

?????文件???????2890??2015-07-16?21:45??三種求解器\DiagImpSym.m

?????文件????????208??2015-07-16?21:08??三種求解器\DiagImpSymCoef.mat

?????文件???????1225??2015-07-16?21:45??三種求解器\EulerCenterForm.m

?????文件????????902??2015-07-17?09:25??三種求解器\Example_1.m

?????文件???????1090??2015-07-16?21:45??三種求解器\R_K_4.m

?????文件???????2432??2015-07-16?21:09??三種求解器\SolveRoots.m

?????目錄??????????0??2015-07-17?09:47??三種求解器

-----------?---------??----------?-----??----

?????????????????8747????????????????????7