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

?屬性????????????大小?????日期????時(shí)間???名稱(chēng)
-----------?---------??----------?-----??----
?????文件????????1833??2018-12-27?00:19??runge4_adaptStep.m
?????文件?????????492??2018-12-27?00:19??test_equationSet.m