在大型语言模型（LLM）驱动的问答和内容生成领域，检索增强生成（Retrieval-Augmented Generation, RAG）已成为一项基础技术。然而，传统的 RAG 依赖于向量相似度搜索，如同在图书馆中寻找内容相似的书籍，却常常忽略了书籍之间、作者之间以及知识点之间错综复杂的关系。为了突破这一局限，Graph RAG 应运而生，它通过引入“知识图谱”，为 RAG 系统装上了能够理解和导航数据关系的“智慧大脑”，实现了从“相关性检索”到“关联性推理”的重大飞跃。

核心概念：从“文档块”到“知识网络”

标准 RAG 的工作流程通常是：

1. 将大量文档切分成独立的文本块（chunks）。

2. 将这些文本块向量化后存入向量数据库。

3. 当用户提问时，将问题向量化并在数据库中进行相似度搜索，找出最相关的几个文本块。

4. 将这些文本块作为上下文（Context）连同问题一起提供给 LLM，生成最终答案。

这种方法的缺陷在于，它处理的是孤立的信息片段。如果一个答案需要整合来自多个不同但相互关联的文档中的信息，标准 RAG 就可能力不从心。

Graph RAG 则彻底改变了数据的组织和检索方式：

1. 构建知识图谱： 它不再仅仅是切分文档，而是利用 LLM 从非结构化文本中提取出关键的实体（Entities）和它们之间的关系（Relationships）。例如，从“埃隆·马斯克是特斯拉的CEO，该公司总部位于美国”这句话中，可以提取出实体“埃隆·G马斯克”、“特斯拉”、“CEO”、“美国”，以及关系“（马斯克）-是-（CEO）-属于-（特斯拉）”和“（特斯拉）-总部位于-（美国）”。

2. 图谱化存储： 这些实体作为“节点（Nodes）”，关系作为“边（Edges）”，共同构建成一个庞大而精细的知识网络，并存储在专门的图数据库中。

3. 基于关系的检索： 当用户提问时，Graph RAG 不仅进行语义搜索，更重要的是在图谱上进行遍历和路径查找。它能识别问题中的核心实体，然后在图谱中探索这些实体的邻近节点和多层关系，从而收集到一个结构化、高度关联且富有深度的上下文。

4. 生成精准答案： 最后，将这个包含丰富实体和关系信息的子图（sub-graph）作为上下文提供给 LLM，使其能够基于清晰的逻辑链条进行推理，生成更精准、更具解释性的答案。

简单来说，标准 RAG 问的是“哪些文档与我的问题最像？”，而 Graph RAG 问的是“与我问题相关的实体有哪些？它们之间存在什么样的联系？这些联系如何共同解答我的问题？”

用到的关键技术

Graph RAG 的实现依赖于一个强大的技术栈，它融合了自然语言处理、图技术和大型模型。

1. 大型语言模型 (LLMs): 在 Graph RAG 中扮演双重角色。

   • 知识提取器： 在构建图谱阶段，使用 LLM 强大的自然语言理解能力，从海量非结构化或半结构化文本中自动识别实体和关系，极大地降低了知识图谱的构建门槛。

   • 推理与生成器： 在查询阶段，LLM 负责理解从图谱中检索到的结构化上下文，并基于这些关系进行推理，最终生成流畅、准确的自然语言答案。

2. 知识图谱 (Knowledge Graphs): 这是 Graph RAG 的核心。它以图的形式组织信息，完美地契合了真实世界中知识相互关联的特性，使得复杂的“多跳问题”（需要跨越多个信息点才能回答的问题）得以解决。

3. 图数据库 (Graph Databases): 专为存储和查询知识图谱而设计的数据库，如 Neo4j, NebulaGraph, Amazon Neptune 等。与传统关系型数据库相比，图数据库能以极高的效率执行关系遍历和路径查找操作，这对于实时响应 Graph RAG 的复杂查询至关重要。

4. 图算法与图机器学习:

   • 社区发现算法 (e.g., Leiden): 在构建图谱后，可以对节点进行聚类，形成不同的“社区”或“主题域”。这使得 Graph RAG 可以实现“全局搜索”，先定位到最相关的社区，再进行细粒度检索。

   • 图嵌入 (Graph Embeddings): 将图中的节点和关系也表示为向量，使得图结构信息可以和文本语义信息在同一个向量空间中进行融合与检索。

主要应用领域

Graph RAG 的优势在于处理需要深度推理、上下文理解和关系挖掘的复杂场景，因此在以下领域展现出巨大潜力：

• 金融风控与合规： 在反欺诈场景中，Graph RAG 可以轻松识别出看似无关的交易、账户和个人之间隐藏的关联，揭示复杂的欺诈网络。在合规审查中，它可以将法规条款、公司行为和交易数据连接起来，提供可追溯、可解释的合规判断。

• 生物医学与药物研发： 生物医学领域充满了复杂的相互作用，如基因、蛋白质、疾病和药物之间的关系。Graph RAG 能够整合海量论文、临床试验数据和基因数据库，帮助研究人员发现新的药物靶点、预测药物副作用，或为特定患者提供个性化的治疗方案建议。

• 供应链管理： 现代供应链是一个复杂的全球网络。Graph RAG 可以整合供应商、物流、库存、市场需求和地缘政治风险等数据，当出现“某国港口关闭”等事件时，能迅速分析出受影响的上下游企业和产品，并推荐替代方案。

• 企业智能知识管理： 大型企业内部知识分散在不同部门的文档、报告和数据库中。Graph RAG 可以将这些知识整合成一个统一的、相互关联的企业知识图谱，让员工能以自然语言提问，获得跨部门、高精度的答案，例如“去年哪个产品的研发项目同时涉及了上海和硅谷的团队，并且项目负责人是谁？”

• 个人化推荐系统： 通过构建用户、产品、偏好和行为之间的关系图谱，Graph RAG 能够提供超越“购买此商品的人也购买了…”的简单推荐，实现更深层次的、基于兴趣和逻辑关联的个性化内容或产品推荐。

What is LangChain
What is LangGraph

LangChain 与 LangGraph：引领大型语言模型应用的双雄

在人工智能应用的浪潮中，大型语言模型（LLM）无疑是当前最耀眼的明星。然而，要将其强大的能力转化为实际、可靠且功能丰富的应用程序，开发者需要一套强而有力的框架。LangChain 和 LangGraph 正是为此而生的两大核心工具，它们极大地简化了 LLM 应用的开发流程，并为构建复杂的智能代理（Agent）提供了无限可能。

LangChain：串联智慧的应用框架

核心概念：

LangChain 是一个开源的应用程序开发框架，其核心理念在于“链（Chain）”的概念。它旨在将大型语言模型与外部资源（如数据库、API、文件等）以及其他计算逻辑无缝地“链接”起来。通过模块化的组件，开发者可以轻松地构建、组合和扩展 LLM 应用，从而实现超越简单文本生成的功能。

LangChain 的主要构成部分包括：

• 模型（Models）： 整合了各类大型语言模型和嵌入模型的标准接口，开发者可以轻松切换和使用不同的模型。

• 提示（Prompts）： 提供强大的提示模板管理功能，让开发者能够动态地生成和优化给予模型的指令。

• 索引（Indexes）： 用于结构化和检索外部数据，是实现如“检索增强生成（Retrieval-Augmented Generation, RAG）”等应用的关键，让模型能够基于特定知识库回答问题。

• 链（Chains）： 这是 LangChain 的核心，用于将多个组件（包括其他链）按顺序组合起来，形成一个完整的应用逻辑。

• 代理（Agents）： 赋予 LLM 推理和行动的能力。代理程序可以使用一系列工具（Tools），并根据用户输入自主决定调用哪个工具来完成任务。

• 记忆（Memory）： 让应用程序能够记住先前的对话内容，从而进行有上下文的、连贯的交流。

应用样例：

LangChain 的应用极为广泛，几乎涵盖了所有需要利用 LLM 进行自然语言处理的场景：

• 智能问答机器人： 结合企业内部文件或特定领域的知识库，打造能够精准回答专业问题的问答系统。例如，一个能够读取公司所有技术文档，并回答工程师技术问题的聊天机器人。

• 文件分析与摘要： 快速读取长篇报告、法律文件或学术论文，并生成简洁的摘要或提取关键信息。

• 数据分析与洞察： 通过自然语言查询，与结构化（如 SQL 数据库）或非结构化数据进行互动，生成分析报告和可视化图表。

• 自动化代码生成与解释： 根据需求描述自动生成代码片段，或对现有代码进行解释和除错。

• 个性化内容推荐： 根据用户的历史行为和偏好，生成个性化的新闻摘要、产品推荐或旅行计划。

LangGraph：驾驭复杂性的图形化流程引擎

核心概念：

随着应用需求的日益复杂，单纯的线性“链”式结构有时会显得力不从心。特别是在需要处理循环、条件分支和多个智能代理协同工作的场景下，LangGraph 应运而生。

LangGraph 是 LangChain 生态系中的一个扩展，它将 LLM 应用的流程从线性的“链”提升到了更灵活、更强大的“图（Graph）”的层次。在 LangGraph 中，应用的执行流程被定义为一个有向图，其中：

• 节点（Nodes）： 代表工作流程中的一个步骤，可以是一个 LLM 的调用、一个工具的执行或一段自定义的 Python 代码。

• 边（Edges）： 连接不同的节点，定义了信息在图中的流动方向。特别是“条件边（Conditional Edges）”，可以根据节点的输出结果，动态地决定下一步要跳转到哪个节点。

• 状态（State）： 这是 LangGraph 的一个核心概念，它是一个在整个图形执行过程中持续存在的对象，用于存储和传递信息。每个节点都可以读取和更新这个共享的状态，从而实现了比线性链更复杂的记忆和上下文管理。

与 LangChain 的主要区别：

应用样例：

LangGraph 特别擅长于构建需要更强控制力和灵活性的高级应用：

• 多代理协同工作（Multi-agent Collaboration）： 构建一个研究团队，其中一个代理负责上网搜索资料，另一个代理负责分析和整合资料，最后一个代理负责撰写报告。LangGraph 可以清晰地定义它们之间的协作流程、信息传递和审核机制。

• 具有自我修正能力的代理： 当一个代理执行的工具出错或结果不理想时，可以通过条件边将流程导向一个“反思与修正”的节点，让代理重新评估策略并尝试新的方法，形成一个不断试错和优化的循环。

• 人机协作（Human-in-the-loop）： 在关键决策点，可以设计一个节点来暂停流程，并请求人类用户的输入或批准。例如，在一个自动下单的代理中，最终下单前需要由人类确认订单细节。

• 互动式故事生成： 根据用户的选择，故事的情节可以走向完全不同的分支，甚至可以回到之前的某个情节点，产生高度互动和非线性的叙事体验。

关键对比：何时用 LangChain，何时用 LangGraph？

目前的发展与未来趋势

LangChain 和 LangGraph 正处于高速发展和迭代的阶段，整个生态系也日益成熟和完善。

LangChain 的发展：

• 架构模块化： LangChain 逐渐将其核心功能拆分为 langchain-core 和 langchain-community 等多个包。langchain-core 提供了稳定和核心的抽象接口，而 langchain-community 则包含了大量的第三方集成，这样的架构使得核心框架更加稳定，同时也方便社区贡献和维护。

• LangChain 表达式语言（LCEL）： 这是 LangChain 近期最重要的更新之一。LCEL 提供了一种声明式的方式来组合链，使得代码更为简洁、直观，并且原生支持批量处理、并行执行、流式输出等高级功能。

• LangSmith 的成熟： LangSmith 是一个用于监控、调试和评估 LLM 应用的平台。它提供了对应用执行过程的“X光”般的洞察力，让开发者可以清晰地看到每一步的输入、输出、延迟和令牌消耗，极大地提升了开发和维护效率。

• LangServe 的推出： 为了方便将开发好的链或代理部署为生产级的 API，LangChain 推出了 LangServe，让开发者可以一键将应用发布为 REST API。

LangGraph 的崛起与未来：

LangGraph 被视为构建下一代复杂智能代理的关键框架。其发展重点在于：

• 成为代理架构的首选： 官方明确表示，LangGraph 将是未来构建代理架构（Agentic Architectures）的重点发展方向。

• 增强可控性和可靠性： LangGraph 的图形化结构和状态管理机制，为设计能够处理复杂任务且行为可预测的代理提供了坚实的基础。

• 人机协作的无缝整合： 其内置的状态持久化和流程控制能力，使得在代理的执行过程中无缝地加入人类的审核和指导成为可能。

• 可视化与开发工具： 配合 LangGraph Studio 等可视化开发工具的推出，开发者可以更直观地设计、调试和迭代复杂的代理流程。

LangChain RAG 代码示例

在实践中，使用 LangChain 来构建一个基于公司内部文档的问答机器人。整个流程非常清晰：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50

# 伪代码和步骤解释 (面试时口述)

# 0. 安装依赖
# pip install langchain langchain-openai langchain-community faiss-cpu

# 1. 加载数据 (Load)
# 我会用 DocumentLoaders (比如 TextLoader 或 PyPDFLoader) 来加载内部文档。
from langchain_community.document_loaders import TextLoader
loader = TextLoader("internal_faq.txt")
documents = loader.load()

# 2. 分割文本 (Split)
# 为了便于检索，需要将长文档切分成小块。
from langchain.text_splitter import RecursiveCharacterTextSplitter
text_splitter = RecursiveCharacterTextSplitter(chunk_size=500, chunk_overlap=0)
docs = text_splitter.split_documents(documents)

# 3. 向量化与存储 (Store)
# 我会使用 OpenAIEmbeddings 将文本块转换为向量，并存入 FAISS 这样的向量数据库中，以便快速搜索。
from langchain_openai import OpenAIEmbeddings
from langchain_community.vectorstores import FAISS
embeddings = OpenAIEmbeddings()
vectorstore = FAISS.from_documents(docs, embeddings)
retriever = vectorstore.as_retriever()

# 4. 创建检索链 (Chain)
# 这是核心步骤。我使用 LangChain 表达式语言 (LCEL) 来构建一个链：
# 它接收问题，使用 retriever 找到相关文档，然后将问题和文档一起传给 LLM 生成答案。
from langchain_core.prompts import ChatPromptTemplate
from langchain_core.runnables import RunnablePassthrough
from langchain_core.output_parsers import StrOutputParser
from langchain_openai import ChatOpenAI

template = """Answer the question based only on the following context:
{context}

Question: {question}
"""
prompt = ChatPromptTemplate.from_template(template)
model = ChatOpenAI()

retrieval_chain = (
    {"context": retriever, "question": RunnablePassthrough()}
    | prompt
    | model
    | StrOutputParser()
)

# 5. 调用链
# print(retrieval_chain.invoke("How to apply for a vacation?"))

“通过这套流程，我能快速搭建一个功能强大的 RAG 系统，整个逻辑是线性和确定的。”

LangGraph 代理代码示例

构建一个更智能的研究代理，它能上网搜索，如果结果不满意，还能决定再次搜索。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62

# 伪代码和步骤解释 (面试时口述)

# 0. 安装依赖
# pip install langgraph langchain-openai tavily-python

# 1. 定义状态 (State)
# 首先，我定义一个图的共享状态，比如包含问题、搜索结果和最终答案。
from typing import TypedDict, List
class AgentState(TypedDict):
    question: str
    search_results: List[str]
    final_answer: str

# 2. 定义节点 (Nodes)
# 每个节点都是一个 Python 函数，接收 state，返回一个更新后的 state 的部分。
def search_tool(state):
    # 调用搜索工具 (如 Tavily)
    # ...
    return {"search_results": ...}

def llm_evaluator(state):
    # 调用 LLM 判断搜索结果是否足够回答问题
    # 如果足够，生成答案
    # 如果不够，返回一个需要重新搜索的信号
    # ...
    return {"final_answer": ...} # 或者返回一个需要再次搜索的标志

# 3. 定义条件边 (Edges)
# 这是 LangGraph 的精髓。我定义一个函数来决定流程走向。
def should_i_search_again(state):
    # 检查 llm_evaluator 的输出
    if "需要重新搜索" in state:
        return "search_tool" # 如果需要，就跳回搜索节点
    else:
        return "__end__" # 否则结束

# 4. 构建图 (Build Graph)
from langgraph.graph import StateGraph, END
graph_builder = StateGraph(AgentState)

# 添加节点
graph_builder.add_node("search_tool", search_tool)
graph_builder.add_node("evaluator", llm_evaluator)

# 设置入口
graph_builder.set_entry_point("search_tool")

# 添加边
graph_builder.add_edge("search_tool", "evaluator")
# 添加条件边
graph_builder.add_conditional_edges(
    "evaluator",
    should_i_search_again,
    {
        "search_tool": "search_tool",
        "__end__": END
    }
)

# 5. 编译和运行
# app = graph_builder.compile()
# app.invoke({"question": "What is the future of AI?"})

For the figure, the each step corresponds to a processed batch(16).

15k dataset

Resdcn + FPNx2 (↓:32,16,8,4)

18 layer

with warm-up and cosine schedule

1
2
3
4
5
6

python train.py --log_name resdcn18_all \
                --data_dir ~/cyl/Data/PSR_final \
                --arch resdcn_18 \
                --lr 1e-5 \
                --batch_size 36 \
                --num_epochs 90 --num_workers 0 --log_interval 5

1
2
3
4
5
6

python train.py --log_name resdcn50_all \
                --data_dir ~/cyl/Data/PSR_final \
                --arch resdcn_50 \
                --lr 1e-5 \
                --batch_size 24 \
                --num_epochs 90 --num_workers 0 --log_interval 5

1
2
3
4
5
6
7

python train.py --log_name resdcn_all \
                --data_dir ~/cyl/Data/PSR_final \
                --arch resdcn_101 \
                --lr 5e-5 \
                --lr_step 40,70 \
                --batch_size 16 \
                --num_epochs 90 --num_workers 0 --log_interval 5

1
2
3
4
5
6
7

python train.py --log_name resdcn_all \
                --data_dir ~/cyl/Data/PSR_final \
                --arch resdcn_101 \
                --lr 5e-5 \
                --lr_step 40,70 \
                --batch_size 16 \
                --num_epochs 90 --num_workers 0 --log_interval 5

HRNet (↓:4,4,4,4)

Use PSR dataset after remove bg samples (11729 samples)

1
2
3
4
5
6
7

python train.py --log_name hrnet_all \
                --data_dir ~/cyl/Data/PSR_final \
                --arch hrnet \
                --lr 5e-5 \
                --lr_step 40,70 \
                --batch_size 16 \
                --num_epochs 90 --num_workers 0 --log_interval 5

1
2
3
4
5
6
7

python train.py --log_name hg_all \
                --data_dir ~/cyl/Data/PSR_final \
                --arch hourglass_small \
                --lr 5e-5 \
                --lr_step 40,70 \
                --batch_size 16 \
                --num_epochs 90 --num_workers 0 --log_interval 5

25k dataset

Resdcn+FPNx2 (↓:32,16,8,4)

50 layer

Resnet (pretrained) (↓:32,16,8,4)

Project Proposal: 2D Point Cloud Registration using the Iterative Closest Point (ICP) Algorithm

Course: VV570/MATH6003J: Introduction to Engineering Numerical Analysis 1

1. Introduction & Problem Statement

In robotic perception, sensors like LiDAR capture the environment as a “point cloud”—a set of data points in space. To enable tasks like mapping or localization, a robot must align multiple point clouds taken from different positions. This project will address this fundamental robotics problem by implementing the Iterative Closest Point (ICP) algorithm, a numerical optimization method used to find the optimal alignment between two point clouds2.

2. Objective

The primary objective of this project is to develop a working 2D implementation of the ICP algorithm. We will apply numerical analysis techniques to find the rigid transformation (rotation and translation) that minimizes the distance between two point clouds, demonstrating a practical application of concepts from the course to the field of robotics.

3. Methodology

Our approach will be to apply a computational tool to solve this problem3. We will use Python with the NumPy library for numerical computations and Matplotlib for visualization. The core steps of our implementation will be:

• Data Association: For each point in a source cloud, find the closest corresponding point in the target cloud.

• Numerical Optimization: Calculate the optimal rotation and translation that minimizes the mean squared error between the corresponding point pairs. This will be solved using Singular Value Decomposition (SVD).

• Iterative Refinement: Apply the calculated transformation to the source cloud and repeat the process until the alignment error converges below a set threshold.

We will start with the simplification of working in 2D and assume a reasonable initial alignment between the point clouds4.

4. Expected Outcomes & Justification of Results

By the end of this three-week project5, we will deliver:

• A Python implementation of the 2D ICP algorithm.

• Visualizations demonstrating the algorithm’s effectiveness by showing the point clouds before and after alignment.

• A plot of the alignment error over iterations, which will serve to justify our results by showing the convergence of the numerical method6.

Homework 5 Report

Zhichen Tang
124370910020

Exercise 1

In this first part, I looked at the logistic model for population growth. The main goal was to solve the equation $u’(t) = C u(t) (1 - u(t)/B)$ using a couple of different numerical methods. The parameters were set to $u_0=100$, $B=1000$, and $C=2/15$.

a)

I take the current value and add the slope at that point multiplied by the time step to get the next value.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

% Completed feuler.m
function [t, u, dt] = feuler( f, I, u0, N )
    dt = (I(2)-I(1))/N;
    t = linspace(I(1),I(2),N+1);
    u = zeros(1,N+1);
    u(1) = u0;
    for n = 1:N
        u(n+1) = u(n) + dt * f(t(n), u(n));
    end
end
````


### b)
First, I calculated a temporary "predicted" value (just like a standard Euler step). Then, I averaged the slope at the current point and the slope at the next point (using the predicted value) to get a more accurate step.

Matlab

% Completed heun.m
function [t, u, dt] = heun( f, I, u0, N )
dt = (I(2)-I(1))/N;
t = linspace(I(1),I(2),N+1);
u = zeros(1,N+1);
u(1) = u0;
for n = 1:N
u_predictor = u(n) + dt * f(t(n), u(n));
u(n+1) = u(n) + (dt/2) * ( f(t(n), u(n)) + f(t(n+1), u_predictor) );
end
end

<div class="post-content"><img src="/2025/07/14/[OBS]课程-Math6003 Hw5/ex1b.png" alt="" title=""></div>
### c) 

After implementing the methods, I ran them with a coarse time step (Deltat=5, corresponding to N=20) and plotted them against the exact solution.

It's pretty clear from the graph that the Heun method does a much better job. Its curve hugs the exact solution line much more tightly than the Forward Euler curve does. This makes sense because it's a higher-order method, so you'd expect it to be more accurate for the same step size.

<div class="post-content"><img src="/2025/07/14/[OBS]课程-Math6003 Hw5/ex1c.png" alt="" title=""></div>

### d)

Next, I repeated the comparison but with a much smaller time step (Deltat=0.05, or N=2000).


Absolutely. With the smaller time step, both methods improved dramatically. The plot shows that the lines for the Euler method, Heun's method, and the exact solution are basically right on top of each other. This really shows how decreasing the step size can significantly reduce the error in numerical solutions. However, the Heun method's absolute error is still lower than the Euler method.

<div class="post-content"><img src="/2025/07/14/[OBS]课程-Math6003 Hw5/ex1d.png" alt="" title=""></div>

---

## Exercise 2

### a)

To solve this, I first had to convert the single 2nd-order equation into a system of two 1st-order equations. I defined a new variable, V(t)=U′(t). This gives the first equation of the system.

Then, I rearranged the original RLC equation to isolate U′′(t):

$$ U''(t) = -\frac{R}{L}U'(t) - \frac{1}{LC}U(t) + \frac{f}{LC} $$

By substituting V for U′ and V′ for U′′, I got the second equation. This allowed me to write the whole system in matrix form, X′=AX+b:

$$\begin{pmatrix} V' \\ U' \end{pmatrix} = \begin{pmatrix} -R/L & -1/(LC) \\ 1 & 0 \end{pmatrix} \begin{pmatrix} V \\ U \end{pmatrix} + \begin{pmatrix} f/(LC) \\ 0 \end{pmatrix} $$
### b)
Stability Analysis of Forward Euler The Forward Euler method is only stable if $|1 + \Delta t \lambda_i| \le 1$ for all eigenvalues ($\lambda_i$) of the matrix A. Using the given component values (L=0.01, C=10, R=0.1), the matrix A becomes: $$ A = \begin{pmatrix} -10 & -10 \\ 1 & 0 \end{pmatrix} $$ I calculated the eigenvalues to be $\lambda = -5 \pm \sqrt{15}$. To ensure stability, the time step $\Delta t$ has to be less than a critical value determined by the eigenvalue with the largest magnitude. The calculation showed that the method is stable only if: $$ \Delta t \le \frac{2}{5 + \sqrt{15}} \approx 0.2254 \text{ s} $$
### c) 
Simulating with Forward Euler I ran simulations for three different time steps. The results were a perfect illustration of the stability condition. 
<div class="post-content"><img src="/2025/07/14/[OBS]课程-Math6003 Hw5/ex2c.png" alt="" title=""></div>
As you can see: - With **N=43** ($\Delta t \approx 0.233$), the time step was too large, and the solution completely blew up. - With **N=46** ($\Delta t \approx 0.217$), the time step was just inside the stable range, and the solution correctly showed a damped wave. - With **N=500** ($\Delta t = 0.02$), the solution was also stable and much smoother. 
### d)
Simulating with Backward Euler Finally, I repeated the same simulations using the Backward Euler method. 
<div class="post-content"><img src="/2025/07/14/[OBS]课程-Math6003 Hw5/ex2d.png" alt="" title=""></div>
The difference is night and day. The Backward Euler method was **stable for all three time steps**, even the large one that made the Forward Euler method fail. This shows that implicit methods like Backward Euler are much more robust for systems like this, as they don't have the same strict stability requirements on the time step. While the accuracy still gets better with smaller steps, you can trust it to give a stable answer no matter what.

使用ssh作为命令行远程工具，启动远程的tensorboardx并且在本地的浏览器中打开。

远程运行：

1

tensorboard --logdir <path> --port 6006

本地运行：

1

ssh -N -f -L localhost:16006:localhost:6006 bohan@10.11.16.146

Ex 1

a

The result of hw4_1_a.m:

Vandermonde matrices are known to be ill-conditioned, with the condition number growing exponentially with $n$.

d

This tridiagonal matrix is well-conditioned. Its condition number grows polynomially with $n$ (like $O(n^2)$). The error $\epsilon_n$ should grow much more slowly than in the Vandermonde case, exhibiting polynomial rather than exponential growth. The residual $r_n$ will be a better, though not perfect, indicator of the true error compared to the previous case.

Comments on the obtained results:

• $\textbf{Error Growth:}$ The relative error $\epsilon_n$ for the tridiagonal system grows polynomially with $n$, which is significantly slower than the exponential growth observed for the Vandermonde matrix. This will be visible as a straight line on the loglog plot.
• $\textbf{Conditioning:}$ The dramatic difference in error behavior is due to the condition number. The Vandermonde matrix is severely ill-conditioned, whereas the tridiagonal matrix is well-conditioned.
• $\textbf{Residual as an Error Indicator:}$ For the well-conditioned tridiagonal matrix, the normalized residual $r_n$ and the relative error $\epsilon_n$ have much closer values. The residual is a far more reliable indicator of the true error in this case than it was for the Vandermonde system.

Ex 2

a

The task is to complete the files $jacobi.m$ and $gauss_seidel.m$. These files implement iterative methods based on a matrix splitting $A = P - (P-A)$, where $P$ is the preconditioner. The iterative update can be expressed as $x^{(k+1)} = x^{(k)} + z^{(k)}$, where $z^{(k)}$ is the solution to $Pz^{(k)} = r^{(k)}$, and $r^{(k)} = b - Ax^{(k)}$ is the residual.

For the $\textbf{Jacobi method}$, the preconditioner P is the diagonal of A. For the $\textbf{Gauss-Seidel method}$, P is the lower triangular part of A. The missing lines implement the calculation of the residual, the update step, and the storing of the residual norm at each iteration.

jacobi.m:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

function [x, iter, res]= jacobi(A,b,x0,nmax,tol)
% JACOBI iterative method
%   [X, Niter, Res] = JACOBI(A,B,X0,NMAX,TOL) attempts to solve the system of linear equations A*X=B
%   for X using the Jacobi method.

% Jacobi preconditioner
P=diag(diag(A)); %

% residual normalization
r0=norm(b);
if r0==0, r0=1; end

% first iteration
x=x0;
r=b-A*x; % The residual b - A*x
res(1)=norm(r); % The norm of the initial residual
iter=0;

while res(end)/r0 > tol && iter < nmax
    z=P\r; % Solve Pz=r
    x=x+z; % Update the solution
    r=b-A*x; % Recompute the residual
    iter=iter+1;
    res(iter+1)=norm(r); % Store the norm of the new residual
end

return

gauss_seidel.m:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

function [x, iter, res]= gauss_seidel(A,b,x0,nmax,tol)
% GAUSS-SEIDEL iterative method
%   [X, Niter, Res] = GAUSS_SEIDEL(A,B,X0,NMAX,TOL) attempts to solve the system of linear equations A*X=B
%   for X using the Gauss-Seidel method.

% Gauss-Seidel preconditioner
P=tril(A); %

% residual normalization
r0=norm(b);
if r0==0, r0=1; end

% first iteration
x=x0;
r=b-A*x; % The residual b - A*x
res(1)=norm(r); % The norm of the initial residual
iter=0;

while res(end)/r0 > tol && iter < nmax
    z=P\r; % Solve Pz=r using forward substitution
    x=x+z; % Update the solution
    r=b-A*x; % Recompute the residual
    iter=iter+1;
    res(iter+1)=norm(r); % Store the norm of the new residual
end

return

b

I solve the system for $n=4$ with the tridiagonal matrix $A$ having $2.5$ on the diagonal and $-1$ on the off-diagonals. The right-hand side is $b=(1.5,0.5,0.5,1.5)^{\top}$. The initial guess is $x^{(0)}=0$, tolerance is $10^{-10}$, and max iterations is $10^8$. I compare the number of iterations and computing time.

The matrix $A$ is strictly diagonally dominant since $|a_{ii}| = 2.5 > \sum_{j \neq i} |a_{ij}|$ for all rows. This property guarantees that both methods will converge.

For this well-conditioned matrix (as shown in part e), the number of iterations grows very slowly with $n$. The condition number remains small and constant. The relative error and the normalized residual are very close in value, indicating that the residual is a good estimator of the error.

d

The analysis from (c) is repeated for the matrix $A$ with 2 on the diagonal and -1 on the off-diagonals , and $b=(1,0,…,0,1)^{\top}$. The exact solution is again $x=(1,…,1)^{\top}$. This matrix is only weakly diagonally dominant, so I expect slower convergence.

The number of iterations will grow much more rapidly with $n$. This is because, as shown in part (e), the condition number of this matrix grows quadratically with $n$, making the problem progressively harder to solve.

e

I use the given formula for the eigenvalues of a symmetric tridiagonal matrix $M$ with $a$ on the diagonal and $b$ on the off-diagonals:
$$ \lambda_{i}(M)=a+2b\cos\left(\frac{\pi i}{n+1}\right), \quad i=1,…,n $$
The condition number is $\kappa_2(A) = |\lambda_{max}|/|\lambda_{min}|$.

Case 1: Matrix from part (c)
Here, $a=2.5$ and $b=-1$. The eigenvalues are $\lambda_i = 2.5 - 2\cos(\frac{\pi i}{n+1})$.

• As $n\to\infty$, $\lambda_{max} \to 2.5 - 2(-1) = 4.5$.
• As $n\to\infty$, $\lambda_{min} \to 2.5 - 2(1) = 0.5$.

The condition number $\kappa(A)$ approaches a constant, $\kappa(A) \to 4.5 / 0.5 = 9$. This is consistent with the numerical results from part (c), explaining the fast and stable convergence.

Case 2: Matrix from part (d)
Here, $a=2$ and $b=-1$. The eigenvalues are $\lambda_i = 2 - 2\cos(\frac{\pi i}{n+1}) = 4\sin^2(\frac{\pi i}{2(n+1)})$.

• As $n\to\infty$, $\lambda_{max} \to 4\sin^2(\pi/2) = 4$.
• For large $n$, $\lambda_{min} = 4\sin^2(\frac{\pi}{2(n+1)}) \approx 4(\frac{\pi}{2(n+1)})^2 = \frac{\pi^2}{(n+1)^2}$.
  The condition number $\kappa(A) \approx \frac{4}{\pi^2/(n+1)^2} = O(n^2)$. This quadratic growth in the condition number is consistent with the numerical results from part (d), explaining the significant increase in iterations with $n$.

Exercise 1

Problem Analysis

The objective is to find the coefficients $(\alpha, \beta, \gamma)$ for the finite difference formula:
$$
D_{h}f(\overline{x})=\frac{\alpha f(\overline{x})+\beta f(\overline{x}-h)+\gamma f(\overline{x}-2h)}{h} \quad
$$
The analysis begins by substituting the Taylor series expansions for $f(\overline{x}-h)$ and $f(\overline{x}-2h)$ around the point $\overline{x}$ into the formula.

The expansions are:
$$
f(\overline{x}-h) = f(\overline{x}) - hf’(\overline{x}) + \frac{h^2}{2}f’’(\overline{x}) - \frac{h^3}{6}f’’’(\overline{x}) + O(h^4)
$$
$$
f(\overline{x}-2h) = f(\overline{x}) - 2hf’(\overline{x}) + 2h^2f’’(\overline{x}) - \frac{4h^3}{3}f’’’(\overline{x}) + O(h^4)
$$
Substituting these into the formula and grouping terms by derivatives of $f(\overline{x})$ results in:
$$
D_{h}f(\overline{x}) = \frac{\alpha + \beta + \gamma}{h} f(\overline{x}) + (-\beta - 2\gamma) f’(\overline{x}) + \left( \frac{\beta}{2} + 2\gamma \right)h f’’(\overline{x}) + O(h^2)
$$
To approximate $f’(\overline{x})$, the coefficient of $f’(\overline{x})$ must be $1$, and the coefficient of $f(\overline{x})$ must be $0$. This yields two necessary equations:

1. $\alpha + \beta + \gamma = 0$
2. $-\beta - 2\gamma = 1$

a) Solution for First-Order Accuracy

For the formula to be of order 1, the two fundamental equations must be satisfied.

From equation (2), $\beta$ can be expressed in terms of $\gamma$:
$$
\beta = -1 - 2\gamma
$$
Substituting this into equation (1):
$$
\alpha + (-1 - 2\gamma) + \gamma = 0 \implies \alpha = 1 + \gamma
$$
Therefore, the family of coefficients $(\alpha, \beta, \gamma)$ that provides at least first-order accuracy is given by:
$$
(\alpha, \beta, \gamma) = (1+\gamma, -1-2\gamma, \gamma), \quad \forall \gamma \in \mathbb{R}
$$

b) Solution for Second-Order Accuracy

To achieve second-order accuracy, the truncation error must be $O(h^2)$. This requires the coefficient of the $h$ term in the error expansion to also be zero.

$$\frac{\beta}{2} + 2\gamma = 0$$

This creates a system of three linear equations to be solved for the unique values of $(\alpha, \beta, \gamma)$:

1. $\alpha + \beta + \gamma = 0$
2. $-\beta - 2\gamma = 1$
3. $\frac{\beta}{2} + 2\gamma = 0$

From equation (3), we find that $\beta = -4\gamma$. Substituting this into equation (2):
$$
-(-4\gamma) - 2\gamma = 1 \implies 2\gamma = 1 \implies \gamma = \frac{1}{2}
$$
With the value for $\gamma$, $\beta$ can be found:
$$
\beta = -4 \left( \frac{1}{2} \right) = -2
$$
Finally, substituting $\beta$ and $\gamma$ into equation (1) yields $\alpha$:
$$
\alpha + (-2) + \frac{1}{2} = 0 \implies \alpha = \frac{3}{2}
$$
The unique values which give a formula of order 2 are:
$$
\left(\alpha, \beta, \gamma\right) = \left(\frac{3}{2}, -2, \frac{1}{2}\right)
$$

Exercise 2

a) Midpoint Formula Calculation

The first task is to compute the integral of the function $f(x)=e^{-x}sin(x)$ on the interval $[a,b]=[0,2]$. The computation uses the composite midpoint formula with $M=10$ sub-intervals.

The resulting numerical approximation of the integral is 0.468592. The exact integral is 0.466630. So the error is 0.001962.

b) Convergence of the Midpoint Formula

This section analyzes the convergence of the composite midpoint formula for the same function and interval. The number of sub-intervals is varied as $M=10^{1},10^{2},10^{3},…10^{5}$.

The absolute error, $|I(f)-Q_{h}^{pm}(f)|$, is computed for each corresponding step size, $h=(b-a)/M$. This error is the difference between the numerical result and the exact value of the integral, $I=(1-e^{-2}(sin(2)+cos(2)))/2$.

A graph of the error versus the step size $h$ is then created using a logarithmic scale for both axes. In such a plot, the slope of the resulting line corresponds to the order of convergence.

The plot displays a straight line with a negative slope. A linear fit to the logarithmic data reveals a slope of approximately 2. This visually confirms that the method exhibits 2nd-order convergence, which aligns with the theoretical error bounds for the composite midpoint rule.

c) Convergence of Trapezoidal and Simpson’s Formulas

The analysis is repeated using the composite trapezoidal formula and the composite Simpson’s formula. The errors for both methods, $|I(f)-Q_{h}^{trap}(f)|$ and $|I(f)-Q_{h}^{simp}(f)|$, are plotted on the same graph in logarithmic scale.

The resulting graph contains two distinct straight lines, one for each method.

• Trapezoidal Rule: The error line for the trapezoidal rule has a slope of approximately 2. This confirms its theoretical 2nd-order convergence.
• Simpson’s Rule: The error line for Simpson’s rule is significantly steeper, with a slope of approximately 4. This demonstrates its theoretical 4th-order convergence.

A comparison of the results shows that for any given step size $h$, the error from Simpson’s rule is substantially smaller than the error from the trapezoidal rule, highlighting the superior accuracy of the higher-order method for smooth functions.

d) Analysis for a Non-Smooth Function

Finally, the convergence analysis is repeated for all three methods on a new function, $f(x)=\sqrt{|x|^{3}}$, over the interval $[a,b]=[-2,2]$. The exact value for this integral is $I=\frac{16}{5}\sqrt{2}$.

Comments on the Obtained Results

The function $f(x)=\sqrt{|x|^{3}}$ is not sufficiently smooth at the point $x=0$ within the integration interval. While the function and its first derivative are continuous, its second derivative, $f’’(x) \propto |x|^{-1/2}$, is unbounded at $x=0$.

The theoretical error estimates that predict 2nd and 4th-order convergence for these methods are based on the assumption that the function’s higher-order derivatives are continuous and bounded. Since this condition is violated, a degradation in the observed convergence rates is expected.

01:

01: