【Runge-Kutta】龍格庫塔不同步長積分到終點不一樣
阿新 • • 發佈:2018-11-19
參考連結:
1、https://blog.csdn.net/xiaokun19870825/article/details/78763739
解析式和微分式
疑問步長不一樣最終結果相同嗎?
對於解析式y=sqrt(1+2x),可以寫成下面常微分形式:
使用RK4(ode45)
下面是自己編寫的matlab程式碼和ode45:
test_fun.m
function dy=test_fun(x,y)
dy = zeros(1,1);%init data
dy(1) = y - 2*x/y;
% dy(1) = 1./y;
runge_kutta1.m
function [x,y]=runge_kutta1(ufunc,y0,h,a,b) % The order of the parameter list is the function name, % initial value vector, step size, time starting point and time % ending point of the differential equation system % (parameter form refers to ode45 function) n=floor((b-a)/h);%Step number x(1)=a;%Starting point of time y(:,1)=y0;%Initial value can be vector, but dimension should be paid attention to. for ii=1:n x(ii+1)=x(ii)+h; k1=ufunc(x(ii),y(:,ii)); k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2); k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2); k4=ufunc(x(ii)+h,y(:,ii)+h*k3); y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6; %Numerical solution based on Runge Kutta method end
執行指令碼run_main.m
[T,F] = ode45(@test_fun1,[0 1],[1]); subplot(131) plot(T,F)%Matlab's own ode45 function results title('Ode45 function effect') [T1,F1]=runge_kutta1(@test_fun1,[1],0.25,0,1);%When testing, change the dimension of test_fun function. Don't forget to change the dimension of the initial value subplot(132) plot(T1,F1)%Self compiled function of Runge Kutta function title('Self compiled Runge Kutta function,step: 0.25') [T2,F2]=runge_kutta1(@test_fun1,[1],0.1,0,1);%When testing, change the dimension of test_fun function. Don't forget to change the dimension of the initial value subplot(133) plot(T2,F2)%Self compiled function of Runge Kutta function title('Self compiled Runge Kutta function,step: 0.01')
從0到1步長不一樣,得到的結果不一樣
ode45:內部自動插值步長0.025,計算出來y(1) = 1.732050816766531
RK4: 步長:0.250,計算出來y(1) = 1.732274190526737
RK4: 步長:0.100, 計算出來y(1) = 1.732056365165567
RK4: 步長:0.025, 計算出來y(1) = 1.732050828604835
RK4: 步長:0.010, 計算出來y(1) = 1.732050808103274 (精度最高)
真實y(1) = 1.732050807568877
結論:針對不含誤差的方程而言:從0積分到1步長越小積分越準確。有誤差的微分方程還未做實驗。