1 %clc, clear
  2 %%求一阶微分方程解析解
  3 % y = dsolve(\'D2y = sqrt(1 + (Dy) ^ 2) / 5 / (1 - x)\', \'y(0) = 0, Dy(0) = 0\', \'x\');
  4 % ezplot(y(2), [0, 0.9999])%explot(S, [x])，显示解析式
  5 % yy = subs(y(2), \'x\', 1)%subs(S, syms, values)
  6 %%求微分方程数值解
  7 % dyy = @(x, yy)[yy(2);sqrt(1 + yy(2) ^ 2) / 5 / (1 - x)];%定义匿名函数
  8 % yy0 = [0, 0]\';%yy(0)和yy(1)的初始值
  9 % [x, yy] = ode45(dyy, [0, 1 - eps], yy0);%求解一阶微分方程的每一个y值;eps:无限接近1
 10 % plot(x, yy(:, 1));%x-y的图象
 11 % yys = yy(end, 1)%x = 1时y的值
 12 %%%
 13 
 14 %%符号解
 15 % y = dsolve(\'x ^ 2 * D2y + x * Dy + (x ^ 2 - 1 / 4) * y\', \'y(pi / 2) = 2, Dy(pi / 2) = - 2 /pi\', \'x\');
 16 % pretty(y)
 17 % ezplot(y)
 18 % hold on
 19 % %数值解
 20 % dyy = @(x, yy) [yy(2); (1 / 4 - x ^ 2) * yy(1) / x ^ 2 - yy(2) / x]; 
 21 % yy0 = [2, -2 / pi]\';
 22 % [x, yy] = ode45(dyy, [pi / 2, 8], yy0);
 23 % plot(x, yy(:, 1), \'*\')
 24 % legend(\'符号解\', \'数值解\')
 25 % %%%
 26 
 27 %%%
 28 %%符号解无解
 29 %%y = dsolve(\'D2y + y * cos(x) = 0\', \'y(0) = 1, Dy(0) = 0\', \'x\')
 30 %ezplot(y)
 31 %%数值解
 32 % yy = @(x) 1 - 1 / gamma(3) * x .^ 2 + 2 / gamma(5) * x .^ 4 - 9 / gamma(7) * x .^ 6 + 55 / gamma(9) * x .^ 8;
 33 % x1 = 0 : 0.1 : 2;
 34 % y1 = yy(x1);
 35 % plot(x1, y1, \'P-\'), hold on %P代表五角星
 36 % dyy = @(x, yy) [yy(2); -yy(1) * cos(x)];
 37 % yy0 = [1, 0]\';
 38 % [x, yy] = ode45(dyy, [0, 2], yy0);
 39 % plot(x, yy(:, 1), \'* -r\')%r代表红色
 40 %legend(\'级数近似解\', \'数值解\', 0)
 41 
 42 %%%
 43 % d = 100; v1= 1; v2 = 2; k = v1 / v2;
 44 % y = @(x) d / 2 * ((x / d) .^ (1 - k) - (x / d) .^ (1 + k));%定义解析解的匿名函数
 45 % ezplot(y, [0, 100])%画解析解的曲线
 46 % dxy = @(t, xy) [-2 * xy(1) / sqrt(xy(1) .^ 2 + xy(2) .^ 2); 1 - 2 * xy(2) / sqrt(xy(1) .^ 2 + xy(2) .^ 2)];%定义微分方程的右端项
 47 % [t, xy] = ode45(dxy, [0, 66.65], [100; 0]);%求数值解，求解的时间区间要逐步试验给出
 48 % solu = [t, xy] %显示数值解
 49 % hold on
 50 % plot(xy(:, 1), xy(:, 2), \'* r\');%画数值解
 51 % legend(\'解析解\', \'数值解\');
 52 % xlabel(\'\'), title(\'\') %不显示x轴标号，不显示标题，标题就是解析式
 53 %%%
 54 
 55 %%隐式微分方程求解
 56 %作为显式微分方程求解
 57 % dx = @(t, x) [-4 * x(1) + x(1) * x(2); x(1) * x(2) - x(2) ^ 2];
 58 % x0 = [2; 1];
 59 % [t, x] = ode45(dx, [0, 10], x0);
 60 % plot(t, x(:, 1), \'-P\');
 61 % hold on;
 62 % plot(t, x(:, 2), \'- *\')
 63 % legend(\'x1的数值解\', \'x2的数值解\');
 64 % [t, x] = ode23(dx, [0, 10], x0);%45比23精确一些
 65 % plot(t, x(:, 1), \'-P\');
 66 % hold on;
 67 % plot(t, x(:, 2), \'- *\')
 68 % legend(\'x1的数值解\', \'x2的数值解\');
 69 %作为隐式微分方程求解
 70 % dxfun = @(t, x, dx) [-dx(1) - 4 * x(1) + x(1) * x(2); -dx(2) + x(1) * x(2) - x(2) ^ 2];
 71 % x0 = [2; 1];
 72 % xp0 = [-6; 1];
 73 % [t, x] = ode15i(dxfun, [0, 10], x0, xp0);
 74 % plot(t, x(:, 1), \'-P\');
 75 % hold on;
 76 % plot(t, x(:, 2), \'-*\');
 77 % legend(\'x1的数值解\', \'x2的数值解\');
 78 %%%
 79 
 80 %%%
 81 %dsolve可以求微分方程组，也可以求二阶的,一般只能求线性（一次方，无互相乘积）常系数（系数是常数，不能形如cos(t)）微分方程,且自由项是某些特殊类型的函数
 82 % [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\')
 83 % y = dsolve(\'y * D2y - Dy ^ 2 = 0\', \'x\')
 84 %%%
 85 
 86 %%求偏导
 87 % syms x y
 88 % f(x, y) = (x + y) / sqrt(x ^ 2 + y ^ 2) * (x ^ 3 - y ^ 3);
 89 % diff(f(x, y), x)
 90 %%黄灯周期与速度关系
 91 % T = 1; a = 10; b = 4.5; mu = 0.2; g = 9.8;
 92 % v0 = [30 50 70] * 1000 / 3600; % 速度单位换算，化成m/s
 93 % y0 = v0 / 2 / mu / g + (a + b) ./ v0 + T
 94 % v = [10 : 75] * 1000 / 3600;% 单位换算
 95 % y = v / 2 / mu / g + (a + b) ./ v + T;
 96 %plot(v, y)
 97 
 98 %%合并符号表达式的同类项
 99 % syms x1 x2 x3;
100 % f1(x1,x2,x3) = x1 + 4 * x2 + 5 * x3;
101 % f2(x1, x2, x3) = x1 - x2 + x3;
102 % f3 = collect(f1 * f2);
103 % diff(f3, x1)
104 
105 %%常微分方程求数值解
106 % F = @(t, y) [y(2); 1000 * (1 - y(1) .^ 2 ) * y(2) - y(1)];
107 % [t, y] = ode15s(F, [0 3000], [2; 0]);
108 % plot(t, y(:, 1), \'o\');
109 % title(\'Solution of van der Po1 Equation, mu = 1000\');
110 % xlabel(\'time t\');
111 % ylabel(\'solution y(:, 1)\');
112 
113 
114 %%常微分方程求解析解
115 % syms x, y;
116 % f = \'D3y - D2y = x\';
117 % dsolve(f, \'y(1) = 8, Dy(1) = 7, D2y(2) = 4\', \'x\')
118 %%常微分方程组求通解和具体解
119 % f1 = \'D2f + 3 * g = sin(x)\';
120 % f2 = \'Dg + Df = cos(x)\';
121 % [f, g] = dsolve(f1, f2, \'x\')
122 % [f, g] = dsolve(f1, f2, \'Df(2) = 0, f(3) = 3, g(5) = 1\', \'x\')
123 
124 %%一阶齐次线性方程组
125 % syms t;
126 % a = [2, 1, 3; 0, 2, -1; 0, 0, 2];
127 % x0 = [1; 2; 1];
128 % x = expm(a * t) * x0
129 
130 %%非齐次线性微分方程组
131 % syms t s;
132 % a = [1, 0, 0; 2, 1, -2; 3, 2, 1];
133 % ft = [0; 0; exp(t) * cos(2*t)];
134 % x0 = [0; 1; 1];
135 % x = expm(a * t) * x0 + int(expm(a * (t - s)) * subs(ft, s), s, 0, t);
136 % %subs函数代值
137 % x = simplify(x)
138 
139 
140 %%int函数求积分int(f, x, x0, xfinal)
141 % syms x
142 % f = 5 / (x - 1) / (x - 2) / (x - 3);
143 % F = int(f, x, 4, 5);
144 % y = eval(F) % eval函数把符号解转化为数值结果
145 
146 %%%边值问题
147 %有对应解的猜测值
148 % yprime = @(x, y) [y(2); (y(1) - 1) * (1 + y(2) ^ 2) ^ (3 / 2)];%f
149 % res = @(ya, yb) [ya(1); yb(1)]; %%% ya(1) = 0, yb(1) = 0边界条件
150 % yinit = @(x) [sqrt(1 - x .^ 2); -x ./ (0.1 + sqrt(1 - x .^ 2))];%边界估计值/初始猜测解
151 % solinit = bvpinit(linspace(-1, 1, 20), yinit);%linspace(start, end, number)
152 % sol = bvp4c(yprime, res, solinit);
153 % fill(sol.x, sol.y(1, :), [0.7, 1, 0.7])
154 % axis([-1, 1, 0, 1])%指定坐标轴范围
155 % xlabel(\'x\', \'FontSize\', 12) %修改标签外观，字体设置为12磅
156 % ylabel(\'h\', \'Rotation\', 0, \'FontSize\', 12)%rotation文本旋转
157 %%%
158 %有参数mu和对应解的猜测值
159 % function sol = skiprun;
160 % solinit = bvpinit(linspace(0, 1, 10), @skipinit, 5);
161 % sol = bvp4c(@skip, @skipbc, solinit);
162 % plot(sol.x, sol.y(1, :), \'-\', sol.x, sol.yp(1, :), \'--\', \'LineWidth\', 2);
163 % xlabel(\'x\', \'FontSize\', 12);
164 % legend(\'y_1\', \'y_2\', 0)
165 % 
166 % function yprime = skip(x, y, mu);
167 % yprime = [y(2); -mu * y(1)];
168 % 
169 % function res = skipbc(ya, yb, mu);
170 % res = [ya(1); ya(2) - 1; yb(1) + yb(2)];
171 % 
172 % function yinit = skipinit(x);
173 % yinit = [sin(x); cos(x)];
174 
175 %%%
176 % function sol = example13;
177 % solinit = bvpinit(linspace(0, 1, 5), @exinit);
178 % sol = bvp4c(@exode, @exbc, solinit); % 函数，边界，猜测初始解
179 % plot(sol.x, sol.y)
180 % 
181 % function yprime = exode(x, y)
182 % yprime = [0.5 * y(1) * (y(3) - y(1)) / y(2)
183 %           -0.5 * (y(3) - y(1))
184 %           (0.9 - 1000 * (y(3) - y(5)) - 0.5 * y(3) * (y(3) - y(1))) / y(4)
185 %           0.5 * (y(3) - y(1))
186 %           100 * (y(3) - y(5))];
187 %       
188 % function res = exbc(ya, yb)
189 % res = [ya(1) - 1; ya(2) - 1; ya(3) - 1; ya(4)+ 10; yb(3) - yb(5)];
190 % 
191 % function yinit = exinit(x)
192 % yinit = [1; 1; -4.5 * x ^ 2 + 8.91 * x + 1; -10; -4.5 * x ^ 2 + 9 * x + 0.91];
193 %%%