【Matlab】抛物方程Crank-Nicolson格式

╰半夏微凉° 2022-12-30 15:51 123阅读 0赞

## 1 题目 ##

(1) 编制用Crank-Nicolson格式求抛物方程数值解的通用程序。  
 \{ ∂ u ∂ t − a ∂ 2 u ∂ x 2 = f ( x , t ) ( 0 < x < 1 , 0 ≤ t ≤   T ) u ( x , 0 ) = φ ( x ) ( 0 ≤ x ≤ 1 ) u ( 0 , t ) = α ( t ) , u ( 1 , t ) = β ( t ) ( 0 < t ≤ T ) \\left\\\{\\begin\{array\}\{ll\} \\frac\{\\partial u\}\{\\partial t\}-a \\frac\{\\partial^\{2\} u\}\{\\partial x^\{2\}\}=f(x, t) & (0<x<1,0 \\leq t \\leq \\mathrm\{~T\}) \\\\ u(x, 0)=\\varphi(x) & (0 \\leq x \\leq 1) \\\\ u(0, t)=\\alpha(t), \\quad u(1, t)=\\beta(t) & (0<t \\leq T) \\end\{array\}\\right. ⎩⎨⎧∂t∂u−a∂x2∂2u=f(x,t)u(x,0)=φ(x)u(0,t)=α(t),u(1,t)=β(t)(0<x<1,0≤t≤ T)(0≤x≤1)(0<t≤T)

(2) 就 a = 1 a=1 a=1, f ( x , t ) = 0 f(x,t)=0 f(x,t)=0, ϕ ( x ) = e x p ( x ) \\phi(x)=exp(x) ϕ(x)=exp(x), α ( t ) = e x p ( t ) \\alpha(t)=exp(t) α(t)=exp(t), β ( t ) = e x p ( t ) \\beta(t)=exp(t) β(t)=exp(t), M = 40 M=40 M=40, N = 40 N=40 N=40,输出点 ( 0.2 , 1.0 ) (0.2,1.0) (0.2,1.0), ( 0.4 , 1.0 ) (0.4,1.0) (0.4,1.0), ( 0.6 , 1.0 ) (0.6,1.0) (0.6,1.0), ( 0.8 , 1.0 ) (0.8,1.0) (0.8,1.0)这四点处 u ( x , t ) u(x,t) u(x,t)的数值解。

(3) 已知所给方程的精确解为 u ( x , t ) = e x p ( x + t ) u(x,t)=exp(x+t) u(x,t)=exp(x\+t)，将步长反复二分，从点 ( 0.2 , 1.0 ) (0.2,1.0) (0.2,1.0), ( 0.4 , 1.0 ) (0.4,1.0) (0.4,1.0), ( 0.6 , 1.0 ) (0.6,1.0) (0.6,1.0), ( 0.8 , 1.0 ) (0.8,1.0) (0.8,1.0)四点处精确解与数值解的误差观察当步长缩小一半时误差以什么规律缩小。

## 2 程序 ##

1 通用程序

function U = CrankNicolson(a,M,N,phi,alpha,beta)
        h = 1 / M;
        tau = 1 / N;
        r = a * tau / h / h;
    
        M1 = zeros(M-1, M-1);
        M2 = zeros(M-1, M-1);
        M3 = zeros(M-1, 1);
        U = zeros(N+1, M-1);
    
        for i=1:M-1
            M1(i,i) = 1 + r;
            M2(i,i) = 1 - r;
            U(1,i) = phi(tau*i);
            if i-1 >= 1
                M1(i, i-1) = - r / 2;
                M2(i, i-1) = r / 2;
            end
            if i+1 <= M-1
                M1(i, i+1) = -r / 2;
                M2(i, i+1) = r / 2;
            end
        end
    
        for k = 1:N
            M3(1) = r /2 * ( alpha(tau*(k-1)) + alpha(tau*k) );
            M3(end) = r /2 * ( beta(tau*(k-1)) + beta(tau*k) );
            temp = inv(M1) *  M2 * U(k,:)' + inv(M1) * M3 ;
            U(k+1,:) = temp';
        end
    end

2 题目二主函数

clc;clear all
    
    %% 初始数据
    a = 1;
    M = 40;
    N = 40; 
    
    syms x t
    phi(x) = exp(x);
    alpha(t) = exp(t);
    beta(t) = exp(1+t);
    
    U = CrankNicolson(a,M,N,phi,alpha,beta);
    
    for j = 0.2:0.2:0.8
        temp = U(end, round(j * M));
        fprintf("输出点（%0.1f,1.0)的数值解为：%f\n",j, temp)
    end

3 题目三主函数

clc;clear all
    
    %% 初始数据
    a = 1;
    
    syms x t
    phi(x) = exp(x);
    alpha(t) = exp(t);
    beta(t) = exp(1+t);
    
    k = 5;
    step = 5;
    
    y = zeros(1,4);
    e = zeros(k,4);
    for j = 1:4
        y(1,j) = exp(0.2 * j + 1);
    end
    
    for i=1:k
        M = step * 2 ^(i-1);
        N = step * 2 ^(i-1); 
        U = CrankNicolson(a,M,N,phi,alpha,beta);
        fprintf("当空间步长为%0.3f时：\n",1 / M)
        for j = 1:4
            temp = U(end, round(j * M / 5));
            e(i,j) = abs(temp - y(1,j));
            fprintf("输出点（%0.1f,1.0) = %f，误差为：%0.8f\n",j, temp,e(i,j))
        end
        fprintf("\n")
        
    end

## 3 运行结果 ##

题目2结果  
![在这里插入图片描述][20201226181011185.png]  
题目3结果  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L1NlZV9TdGFy_size_16_color_FFFFFF_t_70]

[20201226181011185.png]: /images/20221120/a48a4f51d49340e5b244d28231a73f43.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L1NlZV9TdGFy_size_16_color_FFFFFF_t_70]: /images/20221120/6b24bdd12a05467681353121486c44b6.png