rk51test.asv

(有源代码)数值分析作业,本文主要包括两个部分,第一部分是常微分方程(ODE)的三个实验题,第二部分是有关的拓展讨论,包括高阶常微分的求解和边值问题的求解(BVP).文中的算法和算例都是基于Matla

%function rk51test
%runge-kutta problem 51 test
[x1,y1]=rk_wjl('funtest1',0,1,1,0.01,50);
[x2,y2]=rk_wjl('funtest1',0,1,1,0.01,1);
[x3,y3]=rk_wjl('funtest1',0,1,1,0.01,-1);
[x4,y4]=rk_wjl('funtest1',0,1,1,0.01,-50);

y2e=exact51(x2,1);
err2=y2-y2e;
% fprintf('%5.3e\n',err2)
%plot(x2,err2)
y3e=exact51(x3,-1);
err3=y3-y3e;
%  fprintf('%5.3e\n',err3)
y4e=exact51(x4,-50);
err4=y4-y4e;
 fprintf('%5.3e\n',err4)
%  plot(x4,y4e)

figure(1)
plot(x1,y1,'--',x2,y2,'-.',x3,y3,'-..',x4,y4,'-');
legend(['\alpha=50 ';'\alpha=1  ';'\alpha=-1 ';'\alpha=-50']);
title('dy/dt=\alpha*y-\alpha*x+1,\alpha=50,1,-1,-50');
xlabel('x')
ylabel('y')

figure(2)
plot(x1,y1,'--',x2,y2,'-.',x3,y3,'-..',x4,y4,'-');
legend(['\alpha=50 ';'\alpha=1  ';'\alpha=-1 ';'\alpha=-50']);
title('dy/dt=\alpha*y-\alpha*x+1,\alpha=50,1,-1,-50');
ylim([-1,30]);
xlabel('x')
ylabel('y')

%to show the error
y1e=exact51(x1,50);
figure(3)
plot(x1,y1,x1,y1e,'*',x1,y1-y1e)
%ylim([-100,100])
title('The exact solution,\alpha=50,and the error');
legend('\alpha=50','exact solution','error');
xlabel('x')
ylabel('y')

figure(4)
plot(x1,y1,x1,y1e,'.',x1,y1-y1e)
ylim([-100,100])
legend('\alpha=50','exact solution','error');
title('The exact solution,\alpha=50,and the error,Zoom in');
fprintf('%+25.15f\n',y1-y1e)
xlabel('x')
ylabel('y')

figure(5)
error=y1-y1e;
plot(x1(1:20),error(1:20),'*')
title('Error')
xlabel('x')
ylabel('y')


figure(6)
errorbar(x1,y1e,error,':')
title('errorbar')
xlabel('x')
ylabel('y')

[x5,y5]=rk_wjl('funtest1',0,1,1,0.01,-50);
[x6,y6]=rk_wjl('funtest1',0,1,1,0.08,-50);
figure(7)
plot(x5,y5,x6,y6)
title('\alpha=-50,step with 0.01 and 0.08')
ylim([-1,50])
legend('step 0.01','step 0.08')
xlabel('x')
ylabel('y')

index=1:8:80;
%y5(index)
%y6
%fprintf('%25.15f,%25.15f\n',y5(index),y6)
fprintf('%25.15f\n',y5(index))
fprintf('%25.15f\n',y6(1:10))
%fprintf('%25.15f,%25.15f\n',y5(index),y6(1:10))
yy5=exact51(x5,-50);
fprintf('%25.15f\n',yy5(index))
fprintf('%25.15f\n',x5(index))