解线性方程组——最速下降法及图形化表示 | 北太天元 or matlab

news/2024/7/25 21:57:22/文章来源:https://blog.csdn.net/Math_Boy_W/article/details/139042877

一、思路转变

A为对称正定矩阵, A x = b Ax = b Ax=b

求解向量 x x x这个问题可以转化为一个求 f ( x ) f(x) f(x)极小值点的问题,为什么可以这样:

f ( x ) = 1 2 x T A x − x T b + c f(x) = \frac{1}{2}x^TAx - x^Tb + c f(x)=21xTAxxTb+c

可以发现:

∇ f = g r a d f = A x − b \nabla f = \mathrm{grad}f = Ax - b f=gradf=Axb

A A A的正定性可以保证 f ( x ) f(x) f(x)的驻点一定是极小值点。而 A x − b = 0 Ax - b = 0 Axb=0得到的就是 f ( x ) f(x) f(x)的驻点,即:

f ( x ∗ ) = min ⁡ f ( x ) ⇔ ∇ f ( x ∗ ) = A x ∗ − b = 0 f(x^{*}) = \min f(x) \quad \Leftrightarrow \quad \nabla f(x^{*}) = Ax^{*} - b = 0 f(x)=minf(x)f(x)=Axb=0

把解线性方程组的问题,转化为求函数 f ( x ) f(x) f(x)的极小值点。

二、最速下降法

怎么找到这个极小值点?

已知一个多元函数沿其负梯度方向函数值下降得最快。

一种较为形象的解释:

想象自己在半山腰上,要到山脚处:

  • 首先要找好下降方向:负梯度方向
  • 之后沿着选定方向直走
  • 走到不能再下降为止(也就是选定方向的最低点),停下来,再找新的下降方向
  • 重复上面的过程,便能到达山脚

翻译成数学语言

  • 给定任意初值 x 0 x_0 x0,计算残量 r 0 = b − A x 0 r_0 = b - Ax_0 r0=bAx0

  • 选择 P = r 0 P = r_0 P=r0为前进方向,计算:

    α = ( r 0 , r 0 ) ( A r 0 , r 0 ) , x 1 = x 0 + α r 0 \alpha = \frac{\left(r_0, r_0\right)}{\left(Ar_0, r_0\right)}, \quad x_1 = x_0 + \alpha r_0 α=(Ar0,r0)(r0,r0),x1=x0+αr0

  • 重复上面的过程。

算法如下:

三、北太天元 or matlab实现

最速下降法解线性方程组

function [x,k,r] = Gradient_Descent(A,b,x0,e_tol,N)% 最速下降法 解线性方程组% Input: A, b(列向量), x0(初始值),e_tol: error tolerant, N: 限制迭代次数小于 N 次             % Output: x , k(迭代次数), r%   Version:            1.0%   last modified:      2024/05/19n = length(b); k = 0; R = zeros(1,N); % 记录残差r = b - A*x0;x = zeros(n,N); % 记录每次迭代结果x(:,1) = x0;norm_r = norm(r,2); R(1) = norm_r;while norm_r > e_tol && k < Nalpha = r'*r/(r'*A*r);  % 计算步长x(:,k+2) = x(:,k+1) + alpha * r;r = b - A * x(:,k+2); % 残量norm_r = norm(r,2);R(k+2)=norm_r;k = k+1;endx = x(:,1:k+1); % 返回每次的迭代结果r = R(1:k+1);   % 返回每次的残差if k>Nfprintf('迭代超出最大迭代次数');elsefprintf('迭代次数=%i\n',k);end
end

四、数值算例

下面例子中统一 $ N=100,e_tol=10^{-8},x_0 = 0 $

例1

A x = b Ax=b Ax=b
A = [ 4 1 0 0 1 4 1 0 0 1 4 1 0 0 1 4 ] b = [ 6 25 − 11 15 ] A=\begin{bmatrix} 4 & 1 & 0 & 0 \\ 1 & 4 & 1 & 0 \\ 0 & 1 & 4 & 1 \\ 0 & 0 & 1 & 4 \end{bmatrix}\quad b= \begin{bmatrix} 6 \\ 25 \\ -11 \\ 15 \end{bmatrix} A= 4100141001410014 b= 6251115

用最速下降法求 x x x ;

实现

clc;clear all,format long;
N = 100; e_tol = 1e-8;
A = [4, 1, 0, 0; 1, 4, 1, 0;  0, 1, 4, 1; 0, 0, 1, 4];
b = [6; 25; -11; 15];
x0 = [0; 0; 0; 0];
[x11, k1, r11] = Gradient_Descent(A, b, x0, e_tol, N);
x_exact = gsem_column(A, b);% 作图查看误差变化
n = length(b);
k1 = k1 + 1;% 数值解
figure(1);
plot(1:k1, x11(1,:), '-*r', 1:k1, x11(2,:), '-og', 1:k1, x11(3,:), '-+b', 1:k1, x11(4,:), '-dk');
legend('x_1', 'x_2', 'x_3', 'x_4');
title('每个数值解的变化');% 残差变化
figure(2);
plot(1:k1, r11, '-*r');
legend('残差');
title('残差变化');

运行后得到


通过这个例子可以初步看到方法是可行的.

例2

下面这个例子我将形象展示最速下降法的实现特点

A = [3 1; 1 5];
b = [-1;1];
c = 0;

对应函数: f ( x , y ) = 1 2 ( 3 x 2 + 2 ⋅ 1 ⋅ x y + 5 y 2 ) − ( − x + y ) + 0 f(x,y)=\frac{1}{2}\left(3x^2+2\cdot1\cdot xy+5y^2\right)-(-x+y)+0 f(x,y)=21(3x2+21xy+5y2)(x+y)+0

三维表示一下

clc;clear all;format long;
A = [3  1; 1 5]; 
b = [-1;1];
c = 0; 
N = 100; e_tol = 1e-8; x0 =zeros(length(b),1);%x0 =[-0.1;-0.1]x = linspace(-1,1,100); 
y = linspace(-1,1,100);
% 网格化、方便作图
[x, y] = meshgrid(x,y);
% 定义函数 f(x) = 0.5 * x' * A * x - x'*b + c
% 为了作图方便,如下定义
f=@(x,y) 0.5 * (A(1,1) * x.^2 + 2 * A(1,2) * x .* y + A(2,2) * y.^2) - (b(1) * x + b(2) * y) + c;
z = f(x,y);mesh(x,y,z)
[x11,k1,r11] = Gradient_Descent(A,b,x0,e_tol,N);figure(1)
mesh(x,y,z)
hold on
% 绘制最速下降法的每次迭代点
%scatter3(x11(1, :), x11(2, :), f(x11(1,:),x11(2,:)),'r','filled');
plot3(x11(1, :), x11(2, :), f(x11(1,:),x11(2,:)),'r-o');xlabel('x');
ylabel('y');
zlabel('f(x, y)');
title('函数的三维表示');
hold off;

运行后得到

绘制等高线图

figure(2)
hold on 
contour(x,y,z,200)
plot(x11(1, :), x11(2, :), 'r-o');
title('最速下降法特点');
colorbar;

运行后得到


为了展示更清晰,将 $ x_0 $设为 [-0.2;0] ,可以得到这样的图像

由图形可以看出,最速下降法是如何下降的.

从某一点,选定最快的下降方向,下降到不能再下降为止,再重新找新的最快的下降方向.就这样依次进行下去.

由此可以看出最速下降法的优点是容易理解和实现较为简单.

当然也可以看出它还存在很大的改进空间,在每一次选方向时,明明有着更快更好的方向(三角形任意的第三边都更快).


以上图形均在北太天元软件中绘制,matlab同样可以正常运行。

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.luyixian.cn/news_show_1052883.aspx

如若内容造成侵权/违法违规/事实不符,请联系dt猫网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

tty/pty/console/getty/shell/telnet

tty 终端(termimal)= tty(Teletypewriter, 电传打印机),作用是提供一个命令的输入输出环境,在linux下使用组合键ctrl+alt+T打开的就是终端,可以认为terminal和tty是同义词。 tty泛指所有的终端设置,这些是真实存在的设备。 通过tty命令可以查看当前终端连接的设备。…

Linux 一键部署alfresco 6

alfresco 前言 Alfresco是一个流行的企业级开源内容管理系统和协作平台。它提供了丰富的功能,包括文档管理、记录管理、协作工具、工作流管理、搜索和版本控制等。Alfresco还具有灵活的部署选项,可以作为本地部署的软件或云服务来使用。 该平台可以帮助组织管理和存储各种类…

WPS文件没有保存怎么恢复?5个解决方案轻松恢复!

“我在WPS上编辑了一个文件&#xff0c;但是还没来得及将它保存&#xff0c;我不小心就退出软件了&#xff0c;现在不知道有什么方法可以恢复WPS文件呢&#xff1f;大家可以帮帮我吗” WPS作为一款功能强大且用户友好的软件&#xff0c;给我们的工作带来了很多的便利。但我们在…

适用于Android的最佳数据恢复软件

如果您的 Android 设备崩溃&#xff0c;您需要找到一种方法来取回您的数据。幸运的是&#xff0c;有许多数据恢复程序可以帮助您恢复丢失的文件。有些是免费的&#xff0c;而另一些则需要付费。这是适用于Android设备的最佳数据恢复软件列表。 什么是数据恢复软件&#xff1f; …

紫光展锐前沿探索 | 满足未来6G多差异化应用场景的技术体系思考

在6G架构/系统设计中&#xff0c;紫光展锐提出了未来6G空口“一体多翼”的技术体系概念&#xff0c;即“Big-Lite Multi-RAT”。本文将详细对该技术体系展开介绍。 “一体多翼”技术体系通过 “体”&#xff08;Big RAT&#xff09;和“翼”&#xff08;Lite RAT&#xff09;的…

Java语言ADR药物不良反应系统源码Java+IntelliJ+IDEA+MySQL一款先进的药物警戒系统

Java语言ADR药物不良反应系统源码JavaIntelliJIDEAMySQL一款先进的药物警戒系统源码 ADR药物不良反应监测系统是一个综合性的监测平台&#xff0c;旨在收集、报告、分析和评价药品在使用过程中可能出现的不良反应&#xff0c;以确保药品的安全性和有效性。 以下是对该系统的详细…

C#【进阶】俄罗斯方块

俄罗斯方块 文章目录 Test1_场景切换相关BeginScene.csBegionOrEndScene.csEndScene.csGame.csGameScene.csISceneUpdate.cs Test2_绘制对象基类和枚举信息DrawObject.csIDraw.csPosition.cs Test3_地图相关Map.cs Test4_坐标信息类BlockInfo.cs Test5_板砖工人类BlockWorker.…

04_前端三大件JS

文章目录 JavaScript1.JS的组成部分2.JS引入2.1 直接在head中通过一对script标签定义脚本代码2.2创建JS函数池文件&#xff0c;所有html文件共享调用 3.JS的数据类型和运算符4.分支结构5.循环结构6.JS函数的声明7.JS中自定义对象8.JS_JSON在客户端使用8.1JSON串格式8.2JSON在前…

Python与OpenCV:图像处理与计算机视觉实战指南

前言 OpenCV&#xff08;Open Source Computer Vision Library&#xff09;是一个开源的计算机视觉和机器学习软件库&#xff0c;它包含了数百种计算机视觉算法&#xff0c;包括图像处理、视频分析、物体检测、面部识别等。结合Python语言的强大功能&#xff0c;OpenCV可以用于…

java医院管理系统源码(springboot+vue+mysql)

风定落花生&#xff0c;歌声逐流水&#xff0c;大家好我是风歌&#xff0c;混迹在java圈的辛苦码农。今天要和大家聊的是一款基于springboot的医院管理系统。项目源码以及部署相关请联系风歌&#xff0c;文末附上联系信息 。 项目简介&#xff1a; 医院管理系统的主要使用者分…

【AJAX前端框架】Asynchronous Javascript And Xml

1 传统请求及缺点 传统的请求都有哪些&#xff1f; 直接在浏览器地址栏上输入URL。点击超链接提交form表单使用JS代码发送请求 window.open(url)document.location.href urlwindow.location.href url… 传统请求存在的问题 页面全部刷新导致了用户的体验较差。传统的请求导…

Android15 Beta更新速览

Android15 Beta更新速览 前台服务变更 前台服务使应用保持活动状态&#xff0c;以便它们可以执行关键且对用户可见的操作&#xff0c;通常以牺牲电池寿命为代价。在 Android 15 Beta 2 中&#xff0c;dataSync 和 mediaProcessing 前台服务类型现在具有约 6 小时的超时时间&a…

Python 植物大战僵尸

文章目录 效果图项目结构实现思路源代码 效果图 项目结构 实现思路 下面是代码的实现思路&#xff1a; 导入必要的库和模块&#xff1a;首先&#xff0c;我们导入了Python的os、time库以及pygame库&#xff0c;还有植物大战僵尸游戏中用到的各个植物和僵尸的类。 初始化游戏和…

AI写作工具的革命:AIGC如何提升内容生产效率

AIGC&#xff0c;即人工智能生成内容&#xff0c;是一种新兴的内容生产方式&#xff0c;它利用人工智能技术来自动生成文本、图像、音频、视频等多种形式的内容即进入实际应用层面。 所以AI不再是高深的、让人望尘莫及的算力算法&#xff0c;而是真实地贴近了我们的生活&#…

VectorDBBench在windows的调试

VectorDBBench在windows的调试 VectorDBBench是一款向量数据库基准测试工具&#xff0c;支持milvus、Zilliz Cloud、Elastic Search、Qdrant Cloud、Weaviate Cloud 、 PgVector、PgVectorRS等&#xff0c;可以测试其QPS、时延、recall。 VectorDBBench是一款使用python编写的…

鸿蒙ArkUI-X跨语言调用说明:【平台桥接(@arkui-x.bridge)】

平台桥接(arkui-x.bridge) 简介 平台桥接用于客户端&#xff08;ArkUI&#xff09;和平台&#xff08;Android或iOS&#xff09;之间传递消息&#xff0c;即用于ArkUI与平台双向数据传递、ArkUI侧调用平台的方法、平台调用ArkUI侧的方法。 以Android平台为例&#xff0c;Ark…

OM电商系统asp.net

OM电商系统&#xff0c;可以让顾客全面了解商品的详细信息&#xff0c;消除网上购物的信息不对称问题。通过商品分类来组织众多的商品&#xff0c;方便顾客找到所需要的商品。提供客服顾客互动机制&#xff0c;提高顾客的参与度。通过设计合理的订单处理流程&#xff0c;提高顾…

第一个Flutter3项目

配置flutter国内源 首先&#xff0c;配置flutter的国内源&#xff1a; env:PUB_HOSTED_URL"https://pub.flutter-io.cn"; env:FLUTTER_STORAGE_BASE_URL"https://storage.flutter-io.cn"配置gradle国内源 修改gradle\wrapper\gradle-wrapper.properties…

使用Python操作Jenkins

大家好&#xff0c;Python作为一种简洁、灵活且功能丰富的编程语言&#xff0c;可以与各种API轻松集成&#xff0c;Jenkins的API也不例外。借助于Python中的python-jenkins模块&#xff0c;我们可以轻松地编写脚本来连接到Jenkins服务器&#xff0c;并执行各种操作&#xff0c;…

MySQL触发器实战:自动执行的秘密

欢迎来到我的博客&#xff0c;代码的世界里&#xff0c;每一行都是一个故事 &#x1f38f;&#xff1a;你只管努力&#xff0c;剩下的交给时间 &#x1f3e0; &#xff1a;小破站 MySQL触发器实战&#xff1a;自动执行的秘密 前言触发器的定义和作用触发器的定义和作用触发器的…