h = 0.01; % 步长 N = (b-a)/h; ya = 1; % 初值 figure(1) aa = -50:1:0; % a小于0时变化的步长 Rr = []; for i = 1:length(aa) d = aa(i); syms x y d; d = double(aa(i)); f=inline(subs('d*(y-x)+1',{'d'},{d}),'x','y'); R = Rungkuta4(f,a,b,N,ya); Rr = [Rr R]; end [m n ] = size(Rr); n = n/2; k = 1; for i = 1:n x = Rr(:,k); y = Rr(:,k+1); k = k+2; plot(x',y'); hold on; end （问题2）求解程序（Matlab10.0软件） h = 0.01且a=-40 figure(1) a= 0;b= 1; h = 0.01; % 定步长为0.01 N = (b-a)/h; ya = 1; aa = -40; % 定a的值为-40 Rr = []; for i = 1:length(aa) d = aa(i); syms x y d; d = double(aa(i)); f=inline(subs('d*(y-x)+1',{'d'},{d}),'x','y'); R = Rungkuta4(f,a,b,N,ya); Rr = [Rr R]; end [m n ] = size(Rr); n = n/2; k = 1; for i = 1:n x = Rr(:,k); y = Rr(:,k+1); k = k+2; plot(x',y'); hold on; end （问题3）利用函数dsolve求解 clc clear a =-40; syms x y h = 0.01; % 步长 N = (b-a)/h; ya = 1; % 初值 figure(2) aa = 0:1:50; % a大于0时变化的步长 Rr = []; for i = 1:length(aa) d = aa(i); syms x y d; d = double(aa(i)); f=inline(subs('d*(y-x)+1',{'d'},{d}),'x','y'); R = Rungkuta4(f,a,b,N,ya); Rr = [Rr R]; end [m n ] = size(Rr); n = n/2; k = 1; for i = 1:n x = Rr(:,k); y = Rr(:,k+1); k = k+2; plot(x',y'); hold on; end h = 0.2且a=-40 figure(2) a= 0;b= 1; h = 0.2; % 定步长为0.2 N = (b-a)/h; ya = 1; aa = -40; % 定a的值为-40 Rr = []; for i = 1:length(aa) d = aa(i); syms x y d; d = double(aa(i)); f=inline(subs('d*(y-x)+1',{'d'},{d}),'x','y'); R = Rungkuta4(f,a,b,N,ya); Rr = [Rr R]; end [m n ] = size(Rr); n = n/2; k = 1; for i = 1:n x = Rr(:,k); y = Rr(:,k+1); k = k+2; plot(x',y'); hold on; end
y=dsolve('Dy=-10*(y-x)+1','y(0)=1','x') x = 0:0.01:1; yy = subs(y,x); fun=inline('-50*x+50*y+1'); [x,y]=ode15s(fun,[0:0.01:1],1); ys=(a*x - 1)/a + (a*exp(a*x) + 1)/a; plot(x,y,'r',x,ys,'b') B题求解程序 %%M文件： function [t,u,v]=Euler(fun,a,b,n,u0,v0) t=zeros(1,n+1);u=zeros(1,n+1);v=zeros(1,n+1); h=(b-a)/n; t(1)=a;u(1)=u0;v(1)=v0; for k=1:n t(k+1)=t(k)+h; u(k+1)=u(k)+h/2*((-2000*u(k)+999.75*v(k)+1000.25)+... (-2000*u(k+1)+999.75*v(k+1)+1000.25)); v(k+1)=v(k)+h/2*((u(k)-v(k))+u(k+1)-v(k+1)); end %%主程序 fun=inline('-2000u+999.75v+1000.25','u-v','t','u','v'); a=0;b=1;n=10;u0=0;v0=-2;[t,u,v]=Euler(fun,a,b,n,u0,v0) C题求解程序 function yy = Euler(a,b,f,h,y0) %% 改进的Euler syms x y xx = a:h:b; yy(1) = y0; for i = 1:length(xx)-1 x = xx(i); y = yy(i); z1 = subs(f); x = xx(i+1); y = yy(i)+h*z1; z2 = subs(f); yy(i+1) = yy(i)+(h/2)*(z1+z2); end clc clear %% 主程序（C） syms x y f = 2.*y./x+(x.^2).*exp(x); %函数 a = 1;b = 2; % 上下界 h = [0.05 0.1]; %步长 y0 = 0; % 初值 yy1 = Euler(a,b,f,h(1),y0); xx = a:h(1):b; yy = xx.*xx.*((exp(1)).^xx-exp(1)); w1 = abs(yy-yy1)./yy; w1(1) = []; fprintf('步长h=0.05时，数值解、误差如下') yy1 w1 % 步长h=0.05时 yy2 = Euler(a,b,f,h(2),y0); xx = a:h(2):b; yy = xx.*xx.*((exp(1)).^xx-exp(1));
w2 = abs(yy-yy2)./yy; w2(1)=[]; fprintf('步长h=0.1时，数值解、误差如下') yy2 w2 %步长h=0.1时 D题求解程序 function yy = Euler(a,b,f,h,y0) %% Euler法 syms x y xx = a:h:b; yy(1) = y0; for i = 1:length(xx)-1 x = xx(i); y = yy(i); z1 = subs(f); yy(i+1) = yy(i)+h*z1; end clc clear %% 主程序（D） syms x y f = 2.*y./x+(x.^2).*exp(x); %函数 a = 1;b = 2; % 上下界 h = [0.025 0.1]; %步长 y0 = 0; % 初值 yy1 = Euler(a,b,f,h(1),y0); xx = a:h(1):b; yy = xx.*xx.*((exp(1)).^xx-exp(1)); w1 = abs(yy-yy1)./yy; w1(1) = []; fprintf('步长h=0.025时，数值解、误差如下') yy1 w1 % 步长h=0.05时 yy2 = Euler(a,b,f,h(2),y0); xx = a:h(2):b; yy = xx.*xx.*((exp(1)).^xx-exp(1)); w2 = abs(yy-yy2)./yy; w2(1)=[]; fprintf('步长h=0.1时，数值解、误差如下') yy2 w2 %步长h=0.1时
实验总结（由学生填写）：