1.m

来自「利用龙格库塔法（等步长」· M 代码 · 共 47 行

clear all
a=0;b=10;
for j=1:4
    N=10*2^(j-1)-1
    h=(b-a)/(N+1)
    k=0.5;
    v=0;
    c=1;
    p=1;
    q0=50;
    T0=400;
    Tout=300;
    kv=10;
    l=1;m=3;
    xp=a:h:b;
    n=1;
    for x=a+h:h:l
        f(n)=0;
        n=n+1;
    end
    for x=l+h:h:m
    f(n)=q0*sin((x-l)/(m-l)*pi)*2*h^2;
    n=n+1;
    end
    for x=m+h:h:b-h
        f(n)=0
        n=n+1;
    end
    f(1)=(k/h^2+v/2/h)*T0;
    f(n-1)=(2*kv/h-v*kv/k)*Tout;
    A=zeros(N,N);
    A(1,1)=4*k;
    A(1,2)=-2*k-v*p*c*h;
    for i=1:N-1
        A(i+1,i)=-k/h^2-v/2/h;
        A(i,i)=2*k/h^2;
        A(i,i+1)=-k/h^2+v/2/h;
    end
    A(N,N-1)=-2*k/h^2;
    A(N,N)=k/h^2*(2+2*h*kv/k)-v*kv/k
        u=A\f';
    up=[0,u',0]

    plot(xp,up)
    hold on
end