MATLAB解微分方程
一.dsolve函数——常微分方程求解析解
二.龙格——库塔函数（ode45）常微分方程求数值解
三.bvp4c函数——边值问题
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 %%%
