常微分方程求解在MATLAB,MatLab常微分方程及常微分方程组的求解

最近参加了数学建模，对于老师说的Euler算法的不同步长的精度不一样，编写了一个M函数文件来实现这个精度的比较，把函数附上：

function [x,y]= Euler(varargin)

%这里使用可变输出输入函数的

%varargin{1}为求解常微分方程的表达式

%varargin{2}为求解常微分方程的定解条件

%需要给出的变量有常微分方程的范围a，b(varargin{3},varargin{4})

%n为对这个区间的分割(varargin{5})

%xlt写于7月8日

%取得算法需要的变量，并附上容易理解的含义变量

a = varargin{3};

b = varargin{4};%自变量的范围

n = varargin{5};%区间的分割次数

h = (b - a)/n;%步长

Dy = varargin{1}; %常微分方程的表达式

y0 = varargin{2}; %常微分方程的定解条件表达式

%首先求出所给常微分方程问题的精确解

x1 = zeros(n+1,1);

y1 = zeros(n+1,1);

syms f1; syms x;

f1 = dsolve(Dy,y0,'x');

x1(1) = a;

y1(1) = subs(f1,{x},{x1(1)});

for i = 2:(n+1)

x1(i) =

x1(i-1) + h;

y1(i) =

double(subs(f1,{x},{x1(i)}));

end

%利用Euler方法求解近似数值微分解

x2 = zeros(n+1,1);

y2 = zeros(n+1,1);

syms y;

x2(1) = a;

y2(1) = subs(f1,{x},{a});%获得原方程的初解

for i = 2:(n+1)

x2(i) =

x2(i-1) + h;

y2(i) =

y2(i-1) + h .*

double(subs(Dy(5:end),{x,y},{x2(i-1),y2(i-1)}));%特别记录Matlab中的字符串操作，提取子字符串即A(3：6)...

end

%返回经过Euler算法算出x与y的值

x = x2;

y = y2;

%画图进行误差比较

plot(x1,y1,'r');

hold on;

plot(x2,y2,'b');

特此记录，以后写了新的算法再分享