【3. 科学计算案例】Matlab案例代码解析：数值分析篇 |
您所在的位置：网站首页 › matlab1e-3 › 【3. 科学计算案例】Matlab案例代码解析：数值分析篇 |
|
Matlab案例代码解析3. 科学计算案例部分内容 记录于2021-10-24。 DemoBinarySearch：二分法搜索DemoBisection：二分法求解非线性方程DemoEuler：显式 Euler、隐式 Euler、改进 Euler、Runge-kuttaDemoGaussian：Gauss 消元求解线性方程组DemoGaussJordan：Gauss-Jordan 求解线性方程组DemoGold：黄金分割法求解非线性方程DemoLagrange：拉格朗日插值DemoLorenz：洛伦兹曲线DemoNewton：牛顿法求解非线性方程DemoNewtonDown：牛顿下山法求解非线性方程组DemoOptAdvanceRetreat：进退法确定搜索最小值区间DemoOptBFGS：BFGS 算法求解无约束问题DemoOptBroyden：Broyden 族算法求解无约束问题DemoOptFRCG：FR 非线性共轭梯度法求解无约束问题DemoOptGold：黄金分割法求区间最小值点DemoOptNelderMead：单纯形法求解多维无约束问题DemoOptQuadratic：二次插值法求解极值DemoOptReviseNewton：修正牛顿法求解无约束问题DemoOptTrustRegion：信頼域算法实现数据拟合DemoOptTrustRegionBFGS：BFGS 信赖域方法求解函数极小值DemoOptTrustRegionDogleg：信頼域算法求解函数极小值DemoOptTrustRegionNewton：牛顿型信赖域方法求解无约束问题DemoPLSRegress：PLS 回归DemoSearchCriterion：一维线搜索 Armijo, Goldstein, Wolfe 准则DemoSimpson：Simpson 积分DemoSpline：三次样条插值法DemoSSOR：对超松弛迭代法求解线性方程组基础篇 基础Demo篇 智能算法篇 二分法搜索clear;clc; vec = 1:20; indx1 = BinarySearch(vec, 0) indx2 = BinarySearch(vec, 5) indx3 = BinarySearch(vec, 25) function indx = BinarySearch(vec, key) if key = vec(length(vec)) indx = length(vec); else % binary search low = 1; high = length(vec); while low 0 H = H - (H * (dx * dx') * H) / (dx' * H * dx) + (dg * dg') / (dg' * dx); % Hesse矩阵更新 end k = k + 1; x0 = x; end x = x0; end function f = fun(x) f = (x(1) + 10 * x(2))^2 + 5 * (x(3) - 10 * x(4))^2 + (x(2) - 2 * x(3))^2 + 10 * (x(1) - x(4))^2; end function gf = gfun(x) gf = [2 * (x(1) + 10 * x(2)) + 20 * (x(1) - x(4)) 20 * (x(1) + 10 * x(2)) + 2 * (x(2) - 2 * x(3)) 10 * (x(3) - 10 * x(4)) - 4 * (x(2) - 2 * x(3)) -100 * (x(3) - 10 * x(4)) - 20 * (x(1) - x(4))]; endBroyden 族算法求解无约束问题clear;clc; x = OptBroyden([0 0]', 1e-5, 1000) function x = OptBroyden(x0, err, MaxIter) %{ 函数功能：Broyden族算法求解无约束问题: min f(x)； 输入： x0：初始点； err：梯度误差阈值； MaxIter：最大迭代次数； 输出： x：最小值点； %} if nargin = 0.0) dk = -gk; end end if norm(gk) 0 Bk = Bk - (Bk*(dx*dx')*Bk)/(dx'*Bk*dx) + (dg*dg')/(dg'*dx); % Hesse矩阵更新 end end k = k+1; end end function [d, val] = Trust_q(fk, gk, Bk, deltak) %{ 函数功能: 求解信赖域子问题: min qk(d)=fk+gk'*d+0.5*d'*Bk*d, s.t.||d|| Delta*Delta a = Delta / sqrt((du'*du)); elseif dB'*dB > Delta*Delta a = sqrt((Delta*Delta - du'*du) / ((dB - du)'*(dB - du))) + 1; end if a 0.75 && sqrt(abs(d'*d) - Delta) 1 Q_h2(i) = 1 - press(i) / ss(i - 1); else Q_h2(i) = 1; end if Q_h2(i) tol y = A'*x; % 特征向量方向 y = y/norm(y); % 归一化 x0 = x; % 更新 x = A*y; end s = x'*x; Eigenvalue = s; Eigenvector = y; end一维线搜索：Armijo, Goldstein, Wolfe 准则clear;clc; function [alpha, newxk, fk, newfk] = OptArmijo(xk, dk) beta = 0.5; sigma = 0.2; m = 0; mk = 0; maxm = 20; while m = fun(xk) + c2 * lam * gfun(xk)' * dk % 输出最佳的步长 newxk = xk + lam * dk; fk = fun(xk); newfk = fun(newxk); break; % 搜索步长不满足 Goldstein 准则，继续迭代 else a = lam; lam = 0.5 * (a + b); if b = sigma * gfun(xk)' * dk break; else a = alpha; alpha = 0.5 * (a + b); end else b = alpha; alpha = (alpha + a) / 2; end end newxk = xk + alpha * dk; fk = fun(xk); newfk = fun(newxk); end function f = fun(x) f = 100 * (x(1)^2 - x(2))^2 + (x(1) - 1)^2; end function gf = gfun(x) gf = [400 * x(1) * (x(1)^2 - x(2)) + 2 * (x(1) - 1); -200 * (x(1)^2 - x(2))]; endSimpson积分clear;clc; f = @(x) x ./ (1 + sin(x)); I = Simpson(f, 0, pi, 20) function I = Simpson(f, a, b, n) h = (b - a) / n; x = a + (0 : 2*n) * h/ 2; % n = 1 I = h / 6 * (f(x(1)) + 4 * f(x(2)) + f(x(end))); for i = 2:n I = I + h / 6 * (2 * f(x(2*i - 1)) + 4 * f(x(2*i))); end % challenge % function [I, n] = mySimpson(f, a, b, n) % h = (b - a) / n; % x = a + (0 : 2*n) * h/ 2; % % n = 1 % I = h / 6 * (f(x(1)) + 4 * f(x(2)) + f(x(end))); % for i = 2:n % I = I + h / 6 * (2 * f(x(2*i - 1)) + 4 * f(x(2*i))); % end % % challenge % if abs(I - pi) > 1e-10 % [I, n] = mySimpson(f, a, b, n + 1); % end end三次样条插值法clear;clc; x = linspace(0, 3, 5); y = x.^2; xi = linspace(0, 3, 10); yi1= Spline(x, y, xi); yi2 = spline(x, y, xi); xx = linspace(0, 3, 100); yy = xx.^2; plot(xx, yy, xi, yi1, '*', xi, yi2, 'k*'); legend('理论值', 'Spline', 'spline'); title('三次样条插值法'); function yy = Spline(x, y, xx) %{ 函数功能：三次样条插值法； 输入： x：已知点横坐标； y：已知点纵坐标； xx：插值点； 输出： yy：插值点的函数值； %} % = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = n = length(x); a = y(1 : end - 1); b = zeros(n - 1, 1); d = zeros(n - 1, 1); dx = diff(x); dy = diff(y); A = zeros(n); B = zeros(n, 1); % 自然样条端点条件：端点二阶导数为零 A(1, 1) = 1; A(n, n) = 1; for i = 2 : n - 1 A(i, i - 1) = dx(i - 1); A(i, i) = 2*(dx(i - 1) + dx(i)); A(i, i + 1) = dx(i); B(i) = 3*(dy(i) / dx(i) - dy(i - 1) / dx(i - 1)); end c = A \ B; for i = 1 : n - 1 d(i) = (c(i + 1) - c(i)) / (3 * dx(i)); b(i) = dy(i) / dx(i) - dx(i)*(2*c(i) + c(i + 1)) / 3; end m = length(xx); yy = zeros(m, 1); for i = 1 : m for ii = 1 : n - 1 if xx(i) >= x(ii) && xx(i) |
| CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |