程式已經寫了，不過步長你得自己調，當步長較小時，計算時間會很長

另外，tend是時間的終值，你可以設小一些。因為解析解為10*cos(x)，我設成pi，就是計算半個週期。

x""(t)=-x(t)

引入y1=x,y2=x"，則

y1"=y2

y2"=-x=-y1

初始條件為：

y1(0)=10;

y2(0)=0;

將下面兩行百分號之間的內容，儲存成DiffEulerRk4.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function MaxDiffX=DiffEulerRk4(dt,PlotFlag)

%dt是步長

%PlotFlag是否作圖

if nargin

dt=0.01;

end

if nargin

PlotFlag=0;

end

f=inline("[y(2);-y(1)]","t","y"); %微分方程的右邊項

t0=0; %初始時刻

tend=pi; %計算的點數

tt=t0:dt:tend; %一系列離散的點

N=length(tt); %點數

y0=[10;0];

%%(1)尤拉法

EulerY=y0;

for i=2:N

EulerY(:,i)=EulerY(:,i-1)+dt*f(tt(i-1),EulerY(:,i-1));

end

EulerX=EulerY(1,:); %取x

%%(2)龍格庫塔法

RkY=y0;

for i=2:N

k1=f(tt(i-1), RkY(:,i-1));

k2=f(tt(i-1)+dt/2, RkY(:,i-1)+k1*dt/2);

k3=f(tt(i-1)+dt/2, RkY(:,i-1)+k2*dt/2);

k4=f(tt(i-1)+dt, RkY(:,i-1)+k3*dt);

RkY(:,i)=RkY(:,i-1)+(k1+2*k2+2*k3+k4)*dt/6;

end

RkX=RkY(1,:); %取x

%精確解

syms t

analytic=dsolve("D2x=-x","x(0)=10","Dx(0)=0","t");

rightdata=subs(analytic,t,tt);

if PlotFlag

plot(tt,EulerX,"b-",tt,RkX,"r--",tt,rightdata,"g-.")

legend("Euler","Runge-Kutta","analytic")

end

MaxDiffX=[max(abs(RkX-rightdata)),max(abs(EulerX-rightdata))];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

所有題，都得你自己調步長。

輸入：

DiffEulerRk4(0.01,1) %步長取0.01的計算結果，引數為1代表作圖，自己得修改步長

%%下面是變化

Error=[];

Dt=[5e-4,1e-3,2e-3,5e-3,0.01,0.05,0.1];

for dt=Dt %幾種步長，自行修改

dt %檢視dt，步長小，計算量大

Error_1=DiffEulerRk4(dt); %不作圖

Error=[Error;Error_1]; %儲存尤拉法誤差

end

semilogx(Dt,Error)

legend("Euler","RK4")

xlabel("步長")

ylabel("誤差")

title("與理論值誤差")