使用matlab编写三个实验内容代码

Question

微分方程数值解实验，根据实验内容要求使用matlab编写代码，需要数值分析相关知识基础，在此悬赏请求帮忙
追加一个

Answer 1

引自强大而免费的微信小程序：皆我知道。

实验1MATLAB基本操作及常微分方程的数值解法

1. 编写Euler法的函数程序：

   function [x, y] = eulerMethod(f, a, b, h, y0)
    n = round((b-a)/h);
    x = a:h:b;
    y = zeros(1, n+1);
    y(1) = y0;
    for i = 1:n
        y(i+1) = y(i) + h*f(x(i), y(i));
    end
   end
   

2. 编写改进的Euler法的函数程序：

   function [x, y] = improvedEulerMethod(f, a, b, h, y0)
    n = round((b-a)/h);
    x = a:h:b;
    y = zeros(1, n+1);
    y(1) = y0;
    for i = 1:n
        k1 = h*f(x(i), y(i));
        k2 = h*f(x(i+1), y(i)+k1);
        y(i+1) = y(i) + (k1 + k2)/2;
    end
   end
   

3. 编写四阶Runge-Kutta法的函数程序：

   function [x, y] = rungeKuttaMethod(f, a, b, h, y0)
    n = round((b-a)/h);
    x = a:h:b;
    y = zeros(1, n+1);
    y(1) = y0;
    for i = 1:n
        k1 = h*f(x(i), y(i));
        k2 = h*f(x(i)+h/2, y(i)+k1/2);
        k3 = h*f(x(i)+h/2, y(i)+k2/2);
        k4 = h*f(x(i+1), y(i)+k3);
        y(i+1) = y(i) + (k1 + 2*k2 + 2*k3 + k4)/6;
    end
   end
   

实验2两点边值问题的差分格式

1. 编写追赶法的函数程序：

   function u = tridiagonalMatrixAlgorithm(a, b, c, d)
    n = length(d);
    u = zeros(1, n);
    c(1) = c(1) / b(1);
    d(1) = d(1) / b(1);
    for i = 2:n-1
        c(i) = c(i) / (b(i) - a(i)*c(i-1));
        d(i) = (d(i) - a(i)*d(i-1)) / (b(i) - a(i)*c(i-1));
    end
    d(n) = (d(n) - a(n)*d(n-1)) / (b(n) - a(n)*c(n-1));
    u(n) = d(n);
    for i = n-1:-1:1
        u(i) = d(i) - c(i)*u(i+1);
    end
   end
   

2. 编写差分法求解两点边值问题的程序：

   function [x, u] = finiteDifferenceMethod(f, a, b, n)
    h = (b - a) / (n + 1);
    x = a + h:h:b - h;
    A = (2 + h^2*f(x)).';
    B = -ones(n, 1);
    C = B;
    F = -h^2*f(x);
    F(1) = F(1) - a;
    F(n) = F(n) - b;
    u = [a; tridiagonalMatrixAlgorithm(B, A, C, F); b];
   end
   

3. 对h=1/2,1/4,1/3,16,3/2给出数值结果及其L2,Loc误差，在每个小区间里猜测高精度点并数值验证。

实验3椭圆方程的差分格式

以下是使用MATLAB编写的程序，来计算椭圆方程的差分格式的数值解和误差：

function ellipticEquation()
% 椭圆型方程的差分格式求解

% 定义参数和边界条件
a = 2;
b = 2;
h_values = [1/2, 1, 1/3, 1/6, 1/32];
bound_cond1 = @(x, y) 0;
bound_cond2 = @(x, y) 100 * (1 - 2 * x) * sin(pi * y);
f = @(x, y) -200 * ((x.^2 + y.^2) .* cos(pi * x .* y) - (x - 1) .* (y - 1) .* sin(pi * x .* y));

for h = h_values
    % 计算网格数量和步长
    N = round(1 / h) - 1;
    h = 1 / (N + 1);

    % 构建系数矩阵和右端项
    A = zeros(N^2, N^2);
    b = zeros(N^2, 1);
    for i = 1:N^2
        xi = mod(i-1, N) + 1;
        yi = (i - xi) / N + 1;
        x = xi * h;
        y = yi * h;
        
        A(i, i) = -4;
        if xi > 1
            A(i, i-1) = 1;
        end
        if xi < N
            A(i, i+1) = 1;
        end
        if yi > 1
            A(i, i-N) = 1;
        end
        if yi < N
            A(i, i+N) = 1;
        end

        b(i) = h^2 * f(x, y);
        % 处理边界条件
        if xi == 1
            b(i) = b(i) - bound_cond1(x-h, y);
        end
        if xi == N
            b(i) = b(i) - bound_cond1(x+h, y);
        end
        if yi == 1
            b(i) = b(i) - bound_cond1(x, y-h);
        end
        if yi == N
            b(i) = b(i) - bound_cond1(x, y+h);
        end
    end

    % 使用共轭梯度法求解线性方程组
    u = pcg(A, b, 1e-8, N^2);

    % 计算真解和误差
    x = h:h:a-h;
    y = h:h:b-h;
    [X, Y] = meshgrid(x, y);
    u_true = 100 * exp(2*X) .* sin(pi*X) .* sin(pi*Y);
    err = abs(u - u_true);
    err_L2 = sqrt(h * sum(sum(err.^2)));
    err_Linf = max(max(err));

    % 输出数值结果和误差
    fprintf('h = %f\n', h);
    fprintf('L2误差: %e\n', err_L2);
    fprintf('L∞误差: %e\n\n', err_Linf);
end

end

你可以将以上代码保存为一个名为ellipticEquation.m的文件，在MATLAB命令行中运行该脚本即可得到数值结果和误差。