数值计算求初值问题

求几个函数，看课件上，老师推导下面这个函数，注意f(tn,phin)=dphi(tn)/dt(n)，乘一个M之后，Mdphi(tn)/dt(n)=B(tn)phi(tn)+C(tn), B和C的值可以参考老师课件上给的案例，还有要注意给的B的值，是[B(tn),B(tn+1)],还是[B(tn+1),B(tn)]

function [LHS,RHS] = OST(theta,timestep,M,B,C,sol)
LHS=M-theta*timestep*B(1);%[B(tn+1),B(tn)]
RHS=(M+(1-theta)*timestep*B(2))*sol+timestep*(theta*C(1)+(1-theta)*C(2));
end

根据老师推导的函数方法，同理可以推出其他的函数表达式

function [LHS,RHS] = AB2(timestep,M,B,C,sol)
LHS=M;
RHS=(M+3/2*timestep*B(1))*sol(1)+timestep/2*(3*C(1)-B(2)*sol(2)-C(2));%[B(tn),B(tn-1)],
end

function [LHS,RHS] = AM3(timestep,M,B,C,sol)
LHS=M-5/12*timestep*B(1);
RHS=(M+2/3*timestep*B(2))*sol(1)+timestep/12*(5*C(1)+8*C(2)-B(3)*sol(2)-C(3));
end

function [LHS,RHS] = BDF2(timestep,M,B,C,sol)
LHS=3/2*M-timestep*B;
RHS=2*M*sol(1)-1/2*M*sol(2)+timestep*C;
end

clear;
te=0:0.01:2;t=0:0.1:2;p1(1)=0;p2(1)=0;p3(1)=0;dt=0.1;p=exp(-5*te).*(te.*te-2*te+2)-2*exp(-6*te);M=1;C=t.*t.*exp(-5*t);%te和p是原函数的值

for i=1:length(t)-1%要求出所有的值，到i+1就是所有的值
    p1(i+1)=p1(i)+dt*(t(i)*t(i)*exp(-5*t(i))-6*p1(i));
    p2(i+1)=(p2(i)+dt*t(i+1)*t(i+1)*exp(-5*t(i+1)))/(1+6*dt);%注意这边要得到化简，i+1放一边
    p3(i+1)=(p3(i)+1/2*dt*(t(i+1)*t(i+1)*exp(-5*t(i+1))+t(i)*t(i)*exp(-5*t(i))-6*p3(i)))/(1+1/2*dt*6);%这是三种直接求解初值问题的方法
end
sol1=[p3(2) p3(1)];sol2=[p3(2) p3(1)];sol3=[p3(2) p3(1)];p4=[p3(1) p3(2)];p5=[p3(1) p3(2)];p6=[p3(1) p3(2)];B1=[-6 -6];B2=[-6 -6 -6];B3=[-6];
for j=2:length(t)-1
    [L1,R1]=AB2(dt,1,B1,[C(j),C(j-1)],sol1);
    [L2,R2]=AM3(dt,1,B2,[C(j+1),C(j),C(j-1)],sol2);
    [L3,R3]=BDF2(dt,1,B3,[C(j+1)],sol3);
    p4(j+1)=L1\R1;p5(j+1)=L2\R2;p6(j+1)=L3\R3;
    sol1=[p4(j+1),p4(j)];
    sol2=[p5(j+1),p5(j)];
    sol3=[p6(j+1),p6(j)];
end

plot(t,p1,'LineWidth',2),hold on
plot(t,p2,'LineWidth',2),hold on
plot(t,p3,'LineWidth',2),hold on
plot(t,p4,'LineWidth',2),hold on
plot(t,p5,'LineWidth',2),hold on
plot(t,p6,'LineWidth',2),hold on
plot(te,p,'k--','LineWidth',2),hold on
axis([0 2 0 5*10^-3])
legend('Euler Expl.','Euler Impl.','Trapez','AB2','AM3','BDF2','exakt')