MATLAB解微分方程

一.dsolve函数——常微分方程求解析解

二.龙格——库塔函数（ode45）常微分方程求数值解

三.bvp4c函数——边值问题

%clc, clear
%%求一阶微分方程解析解
% y = dsolve('D2y = sqrt(1 + (Dy) ^ 2) / 5 / (1 - x)', 'y(0) = 0, Dy(0) = 0', 'x');
% ezplot(y(2), [0, 0.9999])%explot(S, [x])，显示解析式
% yy = subs(y(2), 'x', 1)%subs(S, syms, values)
%%求微分方程数值解
% dyy = @(x, yy)[yy(2);sqrt(1 + yy(2) ^ 2) / 5 / (1 - x)];%定义匿名函数
% yy0 = [0, 0]';%yy(0)和yy(1)的初始值
% [x, yy] = ode45(dyy, [0, 1 - eps], yy0);%求解一阶微分方程的每一个y值;eps:无限接近1
% plot(x, yy(:, 1));%x-y的图象
% yys = yy(end, 1)%x = 1时y的值
%%%

%%符号解
% y = dsolve('x ^ 2 * D2y + x * Dy + (x ^ 2 - 1 / 4) * y', 'y(pi / 2) = 2, Dy(pi / 2) = - 2 /pi', 'x');
% pretty(y)
% ezplot(y)
% hold on
% %数值解
% dyy = @(x, yy) [yy(2); (1 / 4 - x ^ 2) * yy(1) / x ^ 2 - yy(2) / x]; 
% yy0 = [2, -2 / pi]';
% [x, yy] = ode45(dyy, [pi / 2, 8], yy0);
% plot(x, yy(:, 1), '*')
% legend('符号解', '数值解')
% %%%

%%%
%%符号解无解
%%y = dsolve('D2y + y * cos(x) = 0', 'y(0) = 1, Dy(0) = 0', 'x')
%ezplot(y)
%%数值解
% yy = @(x) 1 - 1 / gamma(3) * x .^ 2 + 2 / gamma(5) * x .^ 4 - 9 / gamma(7) * x .^ 6 + 55 / gamma(9) * x .^ 8;
% x1 = 0 : 0.1 : 2;
% y1 = yy(x1);
% plot(x1, y1, 'P-'), hold on %P代表五角星
% dyy = @(x, yy) [yy(2); -yy(1) * cos(x)];
% yy0 = [1, 0]';
% [x, yy] = ode45(dyy, [0, 2], yy0);
% plot(x, yy(:, 1), '* -r')%r代表红色
%legend('级数近似解', '数值解', 0)

%%%
% d = 100; v1= 1; v2 = 2; k = v1 / v2;
% y = @(x) d / 2 * ((x / d) .^ (1 - k) - (x / d) .^ (1 + k));%定义解析解的匿名函数
% ezplot(y, [0, 100])%画解析解的曲线
% dxy = @(t, xy) [-2 * xy(1) / sqrt(xy(1) .^ 2 + xy(2) .^ 2); 1 - 2 * xy(2) / sqrt(xy(1) .^ 2 + xy(2) .^ 2)];%定义微分方程的右端项
% [t, xy] = ode45(dxy, [0, 66.65], [100; 0]);%求数值解，求解的时间区间要逐步试验给出
% solu = [t, xy] %显示数值解
% hold on
% plot(xy(:, 1), xy(:, 2), '* r');%画数值解
% legend('解析解', '数值解');
% xlabel(''), title('') %不显示x轴标号，不显示标题，标题就是解析式
%%%

%%隐式微分方程求解
%作为显式微分方程求解
% dx = @(t, x) [-4 * x(1) + x(1) * x(2); x(1) * x(2) - x(2) ^ 2];
% x0 = [2; 1];
% [t, x] = ode45(dx, [0, 10], x0);
% plot(t, x(:, 1), '-P');
% hold on;
% plot(t, x(:, 2), '- *')
% legend('x1的数值解', 'x2的数值解');
% [t, x] = ode23(dx, [0, 10], x0);%45比23精确一些
% plot(t, x(:, 1), '-P');
% hold on;
% plot(t, x(:, 2), '- *')
% legend('x1的数值解', 'x2的数值解');
%作为隐式微分方程求解
% dxfun = @(t, x, dx) [-dx(1) - 4 * x(1) + x(1) * x(2); -dx(2) + x(1) * x(2) - x(2) ^ 2];
% x0 = [2; 1];
% xp0 = [-6; 1];
% [t, x] = ode15i(dxfun, [0, 10], x0, xp0);
% plot(t, x(:, 1), '-P');
% hold on;
% plot(t, x(:, 2), '-*');
% legend('x1的数值解', 'x2的数值解');
%%%

%%%
%dsolve可以求微分方程组，也可以求二阶的,一般只能求线性（一次方，无互相乘积）常系数（系数是常数，不能形如cos(t)）微分方程,且自由项是某些特殊类型的函数
% [x, y] = dsolve('Dx + 2 * x - Dy = 10 * cos(t), Dx + Dy + 2 * y = 4 * exp(-2 * t)', 'x(0) = 2, y(0) = 0', 't')
% y = dsolve('y * D2y - Dy ^ 2 = 0', 'x')
%%%

%%求偏导
% syms x y
% f(x, y) = (x + y) / sqrt(x ^ 2 + y ^ 2) * (x ^ 3 - y ^ 3);
% diff(f(x, y), x)
%%黄灯周期与速度关系
% T = 1; a = 10; b = 4.5; mu = 0.2; g = 9.8;
% v0 = [30 50 70] * 1000 / 3600; % 速度单位换算，化成m/s
% y0 = v0 / 2 / mu / g + (a + b) ./ v0 + T
% v = [10 : 75] * 1000 / 3600;% 单位换算
% y = v / 2 / mu / g + (a + b) ./ v + T;
%plot(v, y)

%%合并符号表达式的同类项
% syms x1 x2 x3;
% f1(x1,x2,x3) = x1 + 4 * x2 + 5 * x3;
% f2(x1, x2, x3) = x1 - x2 + x3;
% f3 = collect(f1 * f2);
% diff(f3, x1)

%%常微分方程求数值解
% F = @(t, y) [y(2); 1000 * (1 - y(1) .^ 2 ) * y(2) - y(1)];
% [t, y] = ode15s(F, [0 3000], [2; 0]);
% plot(t, y(:, 1), 'o');
% title('Solution of van der Po1 Equation, mu = 1000');
% xlabel('time t');
% ylabel('solution y(:, 1)');


%%常微分方程求解析解
% syms x, y;
% f = 'D3y - D2y = x';
% dsolve(f, 'y(1) = 8, Dy(1) = 7, D2y(2) = 4', 'x')
%%常微分方程组求通解和具体解
% f1 = 'D2f + 3 * g = sin(x)';
% f2 = 'Dg + Df = cos(x)';
% [f, g] = dsolve(f1, f2, 'x')
% [f, g] = dsolve(f1, f2, 'Df(2) = 0, f(3) = 3, g(5) = 1', 'x')

%%一阶齐次线性方程组
% syms t;
% a = [2, 1, 3; 0, 2, -1; 0, 0, 2];
% x0 = [1; 2; 1];
% x = expm(a * t) * x0

%%非齐次线性微分方程组
% syms t s;
% a = [1, 0, 0; 2, 1, -2; 3, 2, 1];
% ft = [0; 0; exp(t) * cos(2*t)];
% x0 = [0; 1; 1];
% x = expm(a * t) * x0 + int(expm(a * (t - s)) * subs(ft, s), s, 0, t);
% %subs函数代值
% x = simplify(x)


%%int函数求积分int(f, x, x0, xfinal)
% syms x
% f = 5 / (x - 1) / (x - 2) / (x - 3);
% F = int(f, x, 4, 5);
% y = eval(F) % eval函数把符号解转化为数值结果

%%%边值问题
%有对应解的猜测值
% yprime = @(x, y) [y(2); (y(1) - 1) * (1 + y(2) ^ 2) ^ (3 / 2)];%f
% res = @(ya, yb) [ya(1); yb(1)]; %%% ya(1) = 0, yb(1) = 0边界条件
% yinit = @(x) [sqrt(1 - x .^ 2); -x ./ (0.1 + sqrt(1 - x .^ 2))];%边界估计值/初始猜测解
% solinit = bvpinit(linspace(-1, 1, 20), yinit);%linspace(start, end, number)
% sol = bvp4c(yprime, res, solinit);
% fill(sol.x, sol.y(1, :), [0.7, 1, 0.7])
% axis([-1, 1, 0, 1])%指定坐标轴范围
% xlabel('x', 'FontSize', 12) %修改标签外观，字体设置为12磅
% ylabel('h', 'Rotation', 0, 'FontSize', 12)%rotation文本旋转
%%%
%有参数mu和对应解的猜测值
% function sol = skiprun;
% solinit = bvpinit(linspace(0, 1, 10), @skipinit, 5);
% sol = bvp4c(@skip, @skipbc, solinit);
% plot(sol.x, sol.y(1, :), '-', sol.x, sol.yp(1, :), '--', 'LineWidth', 2);
% xlabel('x', 'FontSize', 12);
% legend('y_1', 'y_2', 0)
% 
% function yprime = skip(x, y, mu);
% yprime = [y(2); -mu * y(1)];
% 
% function res = skipbc(ya, yb, mu);
% res = [ya(1); ya(2) - 1; yb(1) + yb(2)];
% 
% function yinit = skipinit(x);
% yinit = [sin(x); cos(x)];

%%%
% function sol = example13;
% solinit = bvpinit(linspace(0, 1, 5), @exinit);
% sol = bvp4c(@exode, @exbc, solinit); % 函数，边界，猜测初始解
% plot(sol.x, sol.y)
% 
% function yprime = exode(x, y)
% yprime = [0.5 * y(1) * (y(3) - y(1)) / y(2)
%           -0.5 * (y(3) - y(1))
%           (0.9 - 1000 * (y(3) - y(5)) - 0.5 * y(3) * (y(3) - y(1))) / y(4)
%           0.5 * (y(3) - y(1))
%           100 * (y(3) - y(5))];
%       
% function res = exbc(ya, yb)
% res = [ya(1) - 1; ya(2) - 1; ya(3) - 1; ya(4)+ 10; yb(3) - yb(5)];
% 
% function yinit = exinit(x)
% yinit = [1; 1; -4.5 * x ^ 2 + 8.91 * x + 1; -10; -4.5 * x ^ 2 + 9 * x + 0.91];
%%%

• » 草稿箱