clear;clc;close all
h=0.2;
t=0:h:3;
x(1)=1;
%使用Runge-Kutta方法，计算微分方程的数值解
for n=1:(length(t)-1)
K1=-t(n)*x(n);
K2=-(t(n)+h/2)*(x(n)+K1*h/2);
K3=-(t(n)+h/2)*(x(n)+K2*h/2);
K4=-(t(n)+h)*(x(n)+K3*h);
x(n+1)=x(n)+(h/6)*(K1+2*K2+2*K3+K4);
end
plot(t,x,'r')
hold on
xt=dsolve('Dx=-t*x','x(0)=1'); %求出微分方程的精确符号解
ezplot(xt,[0 3])
xlabel('t')
ylabel('x')
title(['h=',num2str(h)])
legend({'4th Runge-Kutta','symbolic solution'})
%根据精确符号解，求出在时间t上，x的精确值
n=1;
for t=0:h:3
xt_value(n)=eval(xt);
n=n+1;
end
error=abs(x-xt_value); %Runge-Kutta方法得到的x，与精确的x之间的误差值
figure
t=0:h:3;
plot(t,error) %绘制出 时间t-误差error 图
xlabel('t')
ylabel('error')
msgbox('MATLAB编程答疑，请加QQ: 1530497909','MATLAB答疑','help')
web http://url.cn/NSFcAs -browser
