参考书籍:常用数值算法及其matlab实现,作者:夏省祥
%第10章 常微分方程初值问题的数值解法
%欧拉方法,例10.1
%书籍:常用数值算法及其matlab实现%第10章 常微分方程初值问题的数值解法%欧拉方法function E = eulerdif(fun, a, b, y0, h)%fun:微分方程的右表达式%a, b 为区间%y0 为初值M = floor(b-a)/h ;%离散点的个数M+1T =zeros(1, M+1); Y =zeros(1, M+1); %行向量T = a:h:b;Y(1) = y0;for i = 1:MY(i+1) = Y(i) +h *feval(fun, T(i), Y(i));endE = [T' Y']; %E是M+1行,2列% function E = eulerdif(fun, a, b, y0, h)% %fun:微分方程的右表达式% %a, b 为区间% %y0 为初值% M = floor(b-a)/h ;%离散点的个数M+1% T =zeros( M+1,1); Y =zeros( M+1,1); %列向量% T = a:h:b;% Y(1) = y0;% for i = 1:M%Y(i+1) = Y(i) +h *feval(fun, T(i), Y(i));% end% E = [T Y]; %E是M+1行,2列
%书籍:常用数值算法及其matlab实现%第10章 常微分方程初值问题的数值解法%欧拉方法,例10.1%function E = eulerdif(fun, a, b, y0, h)clear all;clc;close all;%fun =@(x,y)-3*x*y;fun = inline('-3*x*y');a = 0; b =2;y0 = 1;h = 0.1E = eulerdif(fun, a, b, y0, h);plot(E(:,1),E(:,2) ,'r*-'); hold on;exa10_1 = dsolve('Dy = -3*x*y', 'y(0) = 1', 'x'); %求出解析解ezplot(exa10_1, [0 2]); %画出解析解的图像legend('数值解','解析解' )
运行结果如下: