内容发布更新时间 : 2024/11/20 19:34:08星期一 下面是文章的全部内容请认真阅读。

function [x,y1,y2]=my(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y1(1)=y0; y2(1)=y0;

for n=1:length(x)-1

y1(n+1)=y1(n)+h*feval(dyfun,x(n),y1(n)); k1=feval(dyfun,x(n),y2(n));

k2=feval(dyfun,x(n)+h/2,y2(n)+h/2*k1); k3=feval(dyfun,x(n)+h/2,y2(n)+h/2*k2); k4=feval(dyfun,x(n)+h,y2(n)+h*k3); y2(n+1)=y2(n)+h/6*(k1+2*k3+2*k2+k4); end

x=x';y1=y1';y2=y2';

x3=0:0.02:1;y3=1/3*exp(-50*x)+x.^2; plot(x,y1,'r',x,y2,'g',x3,y3,'b'); legend('Euler','Runge-Kutta','真值解')

>> dyfun=inline('-50*y+50*x.^2+2*x'); [x,y1,y2]=my(dyfun,[0,1],1/3,0.02);[x,y1,y2]

H=0.02时

ans =

0 0.3333 0.3333 0.0200 0 0.1254 0.0400 0.0012 0.0485 0.0600 0.0032 0.0212 0.0800 0.0060 0.0130 0.1000 0.0096 0.0125 0.1200 0.0140 0.1400 0.0192 0.1600 0.0252 0.1800 0.0320 0.2000 0.0396 0.2200 0.0480 0.2400 0.0572 0.2600 0.0672 0.2800 0.0780 0.3000 0.0896 0.3200 0.1020 0.3400 0.1152 0.3600 0.1292 0.3800 0.1440 0.4000 0.1596 0.4200 0.1760 0.4400 0.1932 0.4600 0.2112 0.4800 0.2300 0.5000 0.2496 0.5200 0.2700 0.5400 0.2912 0.5600 0.3132 0.5800 0.3360 0.6000 0.3596 0.6200 0.3840 0.6400 0.4092 0.6600 0.4352 0.6800 0.4620 0.7000 0.4896 0.7200 0.5180 0.7400 0.5472 0.7600 0.5772 0.0153 0.0200 0.0257 0.0325 0.0400 0.0484 0.0576 0.0676 0.0784 0.0900 0.1024 0.1156 0.1296 0.1444 0.1600 0.1764 0.1936 0.2116 0.2304 0.2500 0.2704 0.2916 0.3136 0.3364 0.3600 0.3844 0.4096 0.4356 0.4624 0.4900 0.5184 0.5476 0.5776

0.7800 0.6080 0.6084 0.8000 0.6396 0.6400 0.8200 0.6720 0.6724 0.8400 0.7052 0.7056 0.8600 0.7392 0.7396 0.8800 0.7740 0.7744 0.9000 0.8096 0.8100 0.9200 0.8460 0.8464 0.9400 0.8832 0.8836 0.9600 0.9212 0.9216 0.9800 0.9600 0.9604

1.0000 0.9996 1.0000

H=0.01