%自動變步長四階龍格庫塔法求常微分方程的解，微分方程的格式應該是
%y‘（x）=f（xy）;?初始值為（x0y0）
%返回值為x、y兩個向量
function?[xy]=runge4_adaptStep（fx0xNy0）
????x（1）=x0;
????y（1）=y0;
????h=（xN-x0）/100;%假定一個初始步長，在計算中根據情況，這個值可能會被程序更改
????i=1;
????x_temp=x（1）;
????
????stepHalfed=0;%記錄h是否被減半過。對于某些情況，或許會出現以h為步長時步長太大，但以h/2為步長時步長又太小
????????????????????%這種情況下，為了避免步長在h和h/2之間來回跳躍進入死循環，將采用h/2為步長
????e2=1e-6;
????e1=1e-5;
????while?x（i）????????y1=getY（fx_tempy（i）h）;%以h為步長計算，得到x_temp+h處的近似值
????????y11=getY（fx_tempy（i）h/2）;%以h/2為步長連續計算兩次，得到x_temp+h處的近似值
????????y12=getY（fx_temp+h/2y11h/2）;
????????e=abs（（y1-y12）/y1）;%計算上面兩種途徑求得的x+h處的y值之間的比值
????????if?e>e1%步長太大，減小步長后重新計算一輪
????????????h=h/2;
????????????stepHalfed=1;
????????elseif?e>e2%步長合適，保存得到的值并準備下一點的計算
????????????x（i+1）=x_temp+h;
????????????y（i+1）=y1;
????????????x_temp=x（i+1）;
????????????i=i+1;
????????????stepHalfed=0;
????????else%步長太小
????????????if?stepHalfed==1%此時步長是經過縮小后的，如果步長放大，又會過大，會陷入死循環。
????????????????????????????????%因此接受當前步長，保存當前步結果且準備進入下一點的計算
????????????????x（i+1）=x_temp+h;
????????????????y（i+1）=y1;
????????????????x_temp=x（i+1）;
????????????????i=i+1;
????????????????stepHalfed=0;
????????????else%此時步長不是經過縮小后的，確實是太小，需要放大
????????????????h=2*h;%二話不說，把步長放大，從頭再試試????????????????
????????????end
????????end????????
????end
end

function?[newY]=getY（fxyh）
????k1=f（xy）;
????k2=f（x+0.5*hy+0.5*h*k1）;
????k3=f（x+0.5*hy+0.5*h*k2）;
????k4=f（x+hy+h*k3）;
????newY=y+h/6*（k1+2*k2+2*k3+k4）;
end

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----
?????文件????????1833??2018-12-27?00:19??runge4_adaptStep.m
?????文件?????????492??2018-12-27?00:19??test_equationSet.m