四阶龙格库塔和向前欧拉方法和隐式欧拉法以及改进欧拉法

四阶龙格库塔和向前欧拉 方法 快递客服问题件处理详细方法山木方法pdf计算方法pdf华与华方法下载八字理论方法下载 和隐式欧拉法以及改进欧拉法 实验九 欧拉方法 信息与计算科学金融 崔振威 201002034031 一、 实验目的: 1、掌握欧拉算法 设计 领导形象设计圆作业设计ao工艺污水处理厂设计附属工程施工组织设计清扫机器人结构设计 及程序实现 二、 实验内容: 1、p364-9.2.4、p386-9.5.6 三、 实验要求: 主程序: 欧拉方法(前项): function [x,y]=DEEuler(f,a,b,y0,n); %f:一阶常微分方程的一般 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 达式的右端函数 %a:自变量的取值下限 %b:自变量的取值上限 %y0:函数的初值 %n:积分的步数 if nargin<5,n=50; end h=(b-a)/n; x(1)=a;y(1)=y0; for i=1:n x(i+1)=x(i)+h; y(i+1)=y(i)+h*feval(f,x(i),y(i)); end format short 欧拉方法(后项): function [x,y]=BAEuler(f,a,b,y0,n); %f:一阶常微分方程的一般表达式的右端函数 %a:自变量的取值下限 %b:自变量的取值上限 %y0:函数的初值 %n:积分的步数 if nargin<5,n=50; end h=(b-a)/n; x(1)=a;y(1)=y0; for i=1:n x(i+1)=x(i)+h; y1=y(i)+h*feval(f,x(i),y(i)); y(i+1)=y(i)+h*feval(f,x(i+1),y1); end format short 梯形算法: function [I,step,h2] = CombineTraprl(f,a,b,eps) %f 被积函数 %a,b 积分上下限 %eps 精度 %I 积分结果 %step 积分的子区间数 if(nargin ==3) eps=1.0e-4; end n=1; h=(b-a)/2; I1=0; I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h; while abs(I2-I1)>eps n=n+1; h=(b-a)/n; I1=I2; I2=0; for i=0:n-1 x=a+h*i; x1=x+h; I2=I2+(h/2)*(subs(sym(f),findsym(sym(f)),x)+subs(sym(f),findsym(sym(f)),x1)); end end I=I2; step=n; h2=(b-a)/n; 改进欧拉方法: function [x,y]=MoEuler(f,a,b,y0,n); %f:一阶常微分方程的一般表达式的右端函数 %a:自变量的取值下限 %b:自变量的取值上限 %y0:函数的初值 %n:积分的步数 if nargin<5,n=50; end h=(b-a)/n; x(1)=a;y(1)=y0; for i=1:n x(i+1)=x(i)+h; y1=y(i)+h*feval(f,x(i),y(i)); y2=y(i)+h*feval(f,x(i+1),y1); y(i+1)=(y1+y2)/2; end format short 四阶龙格-库塔法: function y = DELGKT4_lungkuta(f, h,a,b,y0,varvec) %f:一阶常微分方程的一般表达式的右端函数 %h:积分步长 %a:自变量的取值下限 %b:自变量的取值上限 %varvec:常微分方程的变量组 format long; N = (b-a)/h; y = zeros(N+1,1); y(1) = y0; x = a:h:b; var = findsym(f); for i=2:N+1 K1 = Funval(f,varvec,[x(i-1) y(i-1)]);%Funval为程序所需要的函数 K2 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]); K3 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K2*h/2]); K4 = Funval(f,varvec,[x(i-1)+h y(i-1)+h*K3]); y(i) = y(i-1)+h*(K1+2*K2+2*K3+K4)/6; end format short; p364-9.2.4 欧拉方法(前项): '2,t21、 y,t,y,y(0),1,y(t),,e，t,2t，2 解: 执行20步时: 编写函数文件doty.m 如下: function f=doty(x,y); f=x^2-y; 在Matlab 命令窗口输入: >> [x1,y1]=DEEuler('doty',0,2,1,20) 可得到结果: x1 = 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 1.1000 1.2000 1.3000 1.4000 1.5000 1.6000 1.7000 1.8000 1.9000 2.0000 y1 = 1.0000 0.9000 0.8110 0.7339 0.6695 0.6186 0.5817 0.5595 0.5526 0.5613 0.5862 0.6276 0.6858 0.7612 0.8541 0.9647 1.0932 1.2399 1.4049 1.5884 1.7906 在Matlab 命令窗口输入: >> y3=-exp(-x1)+x1.^2-2*x1+2 求得解析解: y3 = 1.0000 0.9052 0.8213 0.7492 0.6897 0.6435 0.6112 0.5934 0.5907 0.6034 0.6321 0.6771 0.7388 0.8175 0.9134 1.0269 1.1581 1.3073 1.4747 1.6604 1.8647 输入: >> plot(x1,y1,'o');hold on >> plot(x1,y3,'*');hold on 可得近似值与解析解的图像比较: 执行40步时: 在Matlab 命令窗口输入: >> [x1,y1]=DEEuler('doty',0,2,1,40) 可得到结果: x1 = 0 0.0500 0.1000 0.1500 0.2000 0.2500 0.3000 0.3500 0.4000 0.4500 0.5000 0.5500 0.6000 0.6500 0.7000 0.7500 0.8000 0.8500 0.9000 0.9500 1.0000 1.0500 1.1000 1.1500 1.2000 1.2500 1.3000 1.3500 1.4000 1.4500 1.5000 1.5500 1.6000 1.6500 1.7000 1.7500 1.8000 1.8500 1.9000 1.9500 2.0000 y1 = 1.0000 0.9500 0.9026 0.8580 0.8162 0.7774 0.7417 0.7091 0.6798 0.6538 0.6312 0.6121 0.5967 0.5848 0.5767 0.5724 0.5719 0.5753 0.5826 0.5940 0.6094 0.6290 0.6526 0.6805 0.7126 0.7490 0.7897 0.8347 0.8841 0.9379 0.9961 1.0588 1.1260 1.1977 1.2739 1.3547 1.4401 1.5301 1.6247 1.7240 1.8279 在Matlab 命令窗口输入: >> y3=-exp(-x1)+x1.^2-2*x1+2 求得解析解: y3 = 1.0000 0.9513 0.9052 0.8618 0.8213 0.7837 0.7492 0.7178 0.6897 0.6649 0.6435 0.6256 0.6112 0.6005 0.5934 0.5901 0.5907 0.5951 0.6034 0.6158 0.6321 0.6526 0.6771 0.7059 0.7388 0.7760 0.8175 0.8633 0.9134 0.9679 1.0269 1.0903 1.1581 1.2305 1.3073 1.3887 1.4747 1.5653 1.6604 1.7602 1.8647 输入: >> plot(x1,y1,'o');hold on >> plot(x1,y3,'*');hold on 可得近似值与解析解的图像比较: 从上面结果可以看出，执行40步时近似解的值要接近于解析解，误差更小，结果更精确。 '225、 y,2ty,y(0),1,y(t),1/(1,t) 解: 执行20步时: 编写函数文件doty.m 如下: function f=doty(x,y); f=2*x*y^2; 在Matlab 命令窗口输入: >> [x1,y1]=DEEuler('doty',0,2,1,20) 可得到结果: x1 = 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 1.1000 1.2000 1.3000 1.4000 1.5000 1.6000 1.7000 1.8000 1.9000 2.0000 y1 = 1.0e+176 * 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 6.4573 在Matlab 命令窗口输入: >> y3=1./(1-x1.^2) 求得解析解: y3 = 1.0e+015 * 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 4.5036 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 输入: >> plot(x1,y1,'o');hold on >> plot(x1,y3,'*');hold on 可得近似值与解析解的图像比较: 执行40步时: 在Matlab 命令窗口输入: >> [x1,y1]=DEEuler('doty',0,2,1,40) 可得到结果: x1 = 0.0000 0.0500 0.1000 0.1500 0.2000 0.2500 0.3000 0.3500 0.4000 0.4500 0.5000 0.5500 0.6000 0.6500 0.7000 0.7500 0.8000 0.8500 0.9000 0.9500 1.0000 1.0500 1.1000 1.1500 1.2000 1.2500 1.3000 1.3500 1.4000 1.4500 1.5000 1.5500 1.6000 1.6500 1.7000 1.7500 1.8000 1.8500 1.9000 1.9500 2.0000 y1 = 1.0e+224 * 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 3.6893 Inf Inf Inf Inf Inf Inf Inf Inf Inf在Matlab 命令窗口输入: >> y3=1./(1-x1.^2) 求得解析解: y3 = 1.0e+015 * 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -2.2518 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 输入: >> plot(x1,y1,'o');hold on >> plot(x1,y3,'*');hold on 可得近似值与解析解的图像比较: 欧拉方法(后项): '2,t21、 y,t,y,y(0),1,y(t),,e，t,2t，2 解: 执行20步时: 编写函数文件doty.m 如下: function f=doty(x,y); f=x^2-y; 在Matlab 命令窗口输入: >> [x1,y1]=BAEuler('doty',0,2,1,20) 可得到结果: x1 = 0.0000 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 1.1000 1.2000 1.3000 1.4000 1.5000 1.6000 1.7000 1.8000 1.9000 2.0000 y1 = 1.0000 0.9110 0.8329 0.7665 0.7127 0.6719 0.6449 0.6323 0.6345 0.6520 0.6852 0.7345 0.8003 0.8829 0.9825 1.0995 1.2341 1.3864 1.5567 1.7452 1.9520 在Matlab 命令窗口输入: >> y3=-exp(-x1)+x1.^2-2*x1+2 求得解析解: y3 = 1.0000 0.9052 0.8213 0.7492 0.6897 0.6435 0.6112 0.5934 0.5907 0.6034 0.6321 0.6771 0.7388 0.8175 0.9134 1.0269 1.1581 1.3073 1.4747 1.6604 1.8647 输入: >> plot(x1,y1,'o');hold on >> plot(x1,y3,'*');hold on 可得近似值与解析解的图像比较: 执行40步时: 在Matlab 命令窗口输入: >> [x1,y1]=BAEuler('doty',0,2,1,40) 可得到结果: x1 = 0 0.0500 0.1000 0.1500 0.2000 0.2500 0.3000 0.3500 0.4000 0.4500 0.5000 0.5500 0.6000 0.6500 0.7000 0.7500 0.8000 0.8500 0.9000 0.9500 1.0000 1.0500 1.1000 1.1500 1.2000 1.2500 1.3000 1.3500 1.4000 1.4500 1.5000 1.5500 1.6000 1.6500 1.7000 1.7500 1.8000 1.8500 1.9000 1.9500 2.0000 y1 = 1.0000 0.9526 0.9079 0.8658 0.8267 0.7904 0.7572 0.7272 0.7003 0.6768 0.6566 0.6399 0.6268 0.6172 0.6114 0.6092 0.6109 0.6164 0.6258 0.6392 0.6566 0.6780 0.7035 0.7332 0.7671 0.8052 0.8475 0.8942 0.9451 1.0005 1.0602 1.1243 1.1929 1.2660 1.3435 1.4256 1.5122 1.6034 1.6992 1.7996 1.9046 在Matlab 命令窗口输入: >> y3=-exp(-x1)+x1.^2-2*x1+2 求得解析解: y3 = 1.0000 0.9513 0.9052 0.8618 0.8213 0.7837 0.7492 0.7178 0.6897 0.6649 0.6435 0.6256 0.6112 0.6005 0.5934 0.5901 0.5907 0.5951 0.6034 0.6158 0.6321 0.6526 0.6771 0.7059 0.7388 0.7760 0.8175 0.8633 0.9134 0.9679 1.0269 1.0903 1.1581 1.2305 1.3073 1.3887 1.4747 1.5653 1.6604 1.7602 1.864 输入: >> plot(x1,y1,'o');hold on >> plot(x1,y3,'*');hold on 可得近似值与解析解的图像比较: 从上面结果可以看出，执行40步时近似解的值要接近于解析解，误差更小，结果更精确。 但是后项欧拉要比前项接近解析解。 '225、 y,2ty,y(0),1,y(t),1/(1,t) 解: 执行20步时: 编写函数文件doty.m 如下: function f=doty(x,y); f=2*x*y^2; 在Matlab 命令窗口输入: >> [x1,y1]=BAEuler('doty',0,2,1,20) 可得到结果: x1 = 0.0000 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 1.1000 1.2000 1.3000 1.4000 1.5000 1.6000 1.7000 1.8000 1.9000 2.0000 y1 = 1.0e+119 * 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 4.0544 Inf Inf Inf Inf Inf Inf Inf 在Matlab 命令窗口输入: >> y3=1./(1-x1.^2) 求得解析解: y3 = 1.0e+015 * 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 4.5036 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 输入: >> plot(x1,y1,'o');hold on >> plot(x1,y3,'*');hold on 可得近似值与解析解的图像比较: 执行40步时: 在Matlab 命令窗口输入: >> [x1,y1]=BAEuler('doty',0,2,1,40) 可得到结果: x1 = 0 0.0500 0.1000 0.1500 0.2000 0.2500 0.3000 0.3500 0.4000 0.4500 0.5000 0.5500 0.6000 0.6500 0.7000 0.7500 0.8000 0.8500 0.9000 0.9500 1.0000 1.0500 1.1000 1.1500 1.2000 1.2500 1.3000 1.3500 1.4000 1.4500 1.5000 1.5500 1.6000 1.6500 1.7000 1.7500 1.8000 1.8500 1.9000 1.9500 2.0000 y1 = 1.0e+195 * 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.0834 Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf 在Matlab 命令窗口输入: >> y3=1./(1-x1.^2) 求得解析解: y3 = 1.0e+015 * 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -2.2518 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 输入: >> plot(x1,y1,'o');hold on >> plot(x1,y3,'*');hold on 可得近似值与解析解的图像比较: 梯形算法: '2,t21、 y,t,y,y(0),1,y(t),,e，t,2t，2 解: 在matlab窗口中输入: >> [I,step,h2] = CombineTraprl('-exp(-x1)+x1^2-2*x1+2',0,2) 可得结果: I = 1.8032 step = 29 h2 = 0.0690 由此可以看出，梯形算法求得的结果比执行20步的欧拉算法要接近解析解，但执行40步的 近似解要比梯形算法精确。 '225、 y,2ty,y(0),1,y(t),1/(1,t) 解: 在matlab窗口中输入: >> [I,step,h2] = CombineTraprl('1/(1-x^2)',0,2) 无法得出结果，程序报 Warning: Divide by zero. 错误 改进欧拉方法: '2,t21、 y,t,y,y(0),1,y(t),,e，t,2t，2 解: 执行20步时: 编写函数文件doty.m 如下: function f=doty(x,y); f=x^2-y; 在Matlab 命令窗口输入: >> [x1,y1]=MoEuler('doty',0,2,1,20) 可得到结果: x1 = 0.0000 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 1.1000 1.2000 1.3000 1.4000 1.5000 1.6000 1.7000 1.8000 1.9000 2.0000 y1 = 1.0000 0.9055 0.8219 0.7501 0.6909 0.6450 0.6130 0.5954 0.5929 0.6059 0.6348 0.6800 0.7418 0.8207 0.9167 1.0304 1.1617 1.3111 1.4786 1.6644 1.8687 在Matlab 命令窗口输入: >> y3=-exp(-x1)+x1.^2-2*x1+2 求得解析解: y3 = 1.0000 0.9052 0.8213 0.7492 0.6897 0.6435 0.6112 0.5934 0.5907 0.6034 0.6321 0.6771 0.7388 0.8175 0.9134 1.0269 1.1581 1.3073 1.4747 1.6604 1.8647 输入: >> plot(x1,y1,'o');hold on >> plot(x1,y3,'*');hold on 可得近似值与解析解的图像比较: 执行40步时: 在Matlab 命令窗口输入: >> [x1,y1]=MoEuler('doty',0,2,1,40) 可得到结果: x1 = 0.0000 0.0500 0.1000 0.1500 0.2000 0.2500 0.3000 0.3500 0.4000 0.4500 0.5000 0.5500 0.6000 0.6500 0.7000 0.7500 0.8000 0.8500 0.9000 0.9500 1.0000 1.0500 1.1000 1.1500 1.2000 1.2500 1.3000 1.3500 1.4000 1.4500 1.5000 1.5500 1.6000 1.6500 1.7000 1.7500 1.8000 1.8500 1.9000 1.9500 2.0000 y1 = 1.0000 0.9513 0.9052 0.8619 0.8214 0.7839 0.7494 0.7181 0.6900 0.6652 0.6438 0.6260 0.6116 0.6009 0.5939 0.5907 0.5912 0.5957 0.6040 0.6164 0.6328 0.6532 0.6778 0.7066 0.7395 0.7768 0.8182 0.8641 0.9142 0.9688 1.0277 1.0911 1.1590 1.2313 1.3082 1.3897 1.4756 1.5662 1.6614 1.7612 1.8657 在Matlab 命令窗口输入: >> y3=-exp(-x1)+x1.^2-2*x1+2 求得解析解: y3 = 1.0000 0.9513 0.9052 0.8618 0.8213 0.7837 0.7492 0.7178 0.6897 0.6649 0.6435 0.6256 0.6112 0.6005 0.5934 0.5901 0.5907 0.5951 0.6034 0.6158 0.6321 0.6526 0.6771 0.7059 0.7388 0.7760 0.8175 0.8633 0.9134 0.9679 1.0269 1.0903 1.1581 1.2305 1.3073 1.3887 1.4747 1.5653 1.6604 1.7602 1.8647 输入: >> plot(x1,y1,'o');hold on >> plot(x1,y3,'*');hold on 可得近似值与解析解的图像比较: 从上面结果可以看出，改进的欧拉法更接近于解析解，并且执行40步时近似解的值要接近于解析解，误差更小，结果更精确。 '225、 y,2ty,y(0),1,y(t),1/(1,t) 解: 执行20步时: 编写函数文件doty.m 如下: function f=doty(x,y); f=2*x*y^2; 在Matlab 命令窗口输入: >> [x1,y1]=MoEuler('doty',0,2,1,20) 可得到结果: x1 = 0.0000 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 1.1000 1.2000 1.3000 1.4000 1.5000 1.6000 1.7000 1.8000 1.9000 2.0000 y1 = 1.0e+113 * 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 3.6181 Inf Inf Inf Inf Inf Inf 在Matlab 命令窗口输入: >> y3=1./(1-x1.^2) 求得解析解: y3 = 1.0e+015 * 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 4.5036 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 输入: >> plot(x1,y1,'o');hold on >> plot(x1,y3,'*');hold on 可得近似值与解析解的图像比较: 执行40步时: 在Matlab 命令窗口输入: >> [x1,y1]=MoEuler('doty',0,2,1,40) 可得到结果: x1 = 0.0000 0.0500 0.1000 0.1500 0.2000 0.2500 0.3000 0.3500 0.4000 0.4500 0.5000 0.5500 0.6000 0.6500 0.7000 0.7500 0.8000 0.8500 0.9000 0.9500 1.0000 1.0500 1.1000 1.1500 1.2000 1.2500 1.3000 1.3500 1.4000 1.4500 1.5000 1.5500 1.6000 1.6500 1.7000 1.7500 1.8000 1.8500 1.9000 1.9500 2.0000 y1 = 1.0e+102 * 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 2.1796 Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf 在Matlab 命令窗口输入: >> y3=1./(1-x1.^2) 求得解析解: y3 = 1.0e+015 * 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -2.2518 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 输入: >> plot(x1,y1,'o');hold on >> plot(x1,y3,'*');hold on 可得近似值与解析解的图像比较: 四阶龙格-库塔法 '2,t21、 y,t,y,y(0),1,y(t),,e，t,2t，2 解: 执行20步时: 在matlab窗口中输入: >> syms x y; >> z=x^2-y >> y= DELGKT4_lungkuta(z,0.1,0,2,1,[x y]) 可以得出下面结果: y = 1.0000 0.9052 0.8213 0.7492 0.6897 0.6435 0.6112 0.5934 0.5907 0.6034 0.6321 0.6771 0.7388 0.8175 0.9134 1.0269 1.1581 1.3073 1.4747 1.6604 1.8647 执行40步时 在matlab窗口中输入: >> syms x y; >> z=x^2-y 可以得出下面结果 >> y= DELGKT4_lungkuta(z,0.05,0,2,1,[x y]) y = 1.0000 0.9513 0.9052 0.8618 0.8213 0.7837 0.7492 0.7178 0.6897 0.6649 0.6435 0.6256 0.6112 0.6005 0.5934 0.5901 0.5907 0.5951 0.6034 0.6158 0.6321 0.6526 0.6771 0.7059 0.7388 0.7760 0.8175 0.8633 0.9134 0.9679 1.0269 1.0903 1.1581 1.2305 1.3073 1.3887 1.4747 1.5653 1.6604 1.7602 1.8647 由上面结果可以看出，无论执行40步还是20步，所求得的近似解相同，但和欧拉法(向 前、向后)相比，该算法更精确，误差更小。 '225、 y,2ty,y(0),1,y(t),1/(1,t) 解: 执行20步时: 在matlab窗口中输入: >> syms x y; >> z=2*x*y^2; >> y= DELGKT4_lungkuta(z,0.1,0,2,1,[x y]) 可以得出下面结果 : y = 1.0e+180 * 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 3.8782 Inf Inf Inf Inf Inf Inf Inf Inf 执行40步时: 在matlab窗口中输入: >> syms x y; >> z=2*x*y^2; >> y= DELGKT4_lungkuta(z,0.05,0,2,1,[x y]) 可以得出以下结果: y = 1.0e+176 * 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 7.4255 Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf 从上面所有算法可以看出，改进的欧拉法要比欧拉法(向前、向后)的好，求得的结果要更 接近解析解，并且执行40步时比执行20步时要误差要小，结果更好。同样四阶龙格-库塔 比欧拉法好，但是不管执行40步还是20步结果无明显差别，因此四阶龙格-库塔法的结果 与执行步数无关。