%书籍：常用数值算法及其matlab实现%第10章 常微分方程初值问题的数值解法%改进欧拉方法function S = impeuler(fun, x0, xn, y0, h)%fun:微分方程的右表达式%x0, xn为区间%y0 为初值M = floor(xn-x0)/h ;%离散点的个数M+1T =zeros(1, M+1); Y =zeros(1, M+1); %行向量T(1) = x0;Y(1) = y0;for i = 1:MK1 = feval(fun, T(i) ,Y(i));K2 = feval(fun, T(i)+h ,Y(i)+ h*K1);Y(i+1) = Y(i) +h/2 *(K1 + K2);T(i+1) = T(i) +h;endS = [T' Y']; %E是M+1行，2列画向量%书籍：常用数值算法及其matlab实现%第10章 常微分方程初值问题的数值解法%欧拉方法，例10.1%画出微分方程的向量场function[] = direction_field(f, a, b, c, d)%画出一阶微分方程y' = f(x,y)的向量场，a<= x <= b, c<= y <= dt1 = a:0.2:b;t2 = c:0.2:d;[x,y] = meshgrid(t1,t2);S = f(x,y);K = inline(vectorize('1/sqrt(1 + S*S)'), 'S');L =K(S)quiver(x, y, L, S.*L, 0.5)axis equal tight主函数%书籍：常用数值算法及其matlab实现%第10章 常微分方程初值问题的数值解法%改进欧拉方法，例10.2%function S = impeuler(fun, x0, xn, y0, h)clear all;clc;close all;%fun =@(x,y)y*cos(x);fun = inline('y*cos(x)');a = 0; b =3;y0 = 1;h = 0.2S = impeuler(fun, a, b, y0, h)plot(S(:,1),S(:,2) ,'r*-'); hold on;exa10_2 = dsolve('Dy = y*cos(x)', 'y(0) = 1', 'x'); %求出解析解ezplot(exa10_2, [0 3]); %画出解析解的图像legend('数值解','解析解' )figure(2) %画出向量场example10_2 = inline(vectorize('y*cos(x)'),'x','y');direction_field(example10_2, 0, 3,0,3); %画出向量场hold on;ezplot(exa10_2, [0 3]); %画出解析解的图像运行结果：
