本学期微分方程数值解课程总结（matlab代码）

Posted 2020-07-21 视野

1. 最简单求解一个微分方程数值解：Euler法

function [x,y]=Euler_method(dufun,span,h,x0,y0)
%EuLer格式,
%求解方程y‘=dufun(x,y);其中x \in[a,b];y0为初始值；n为自变量的离散个数；y为求解结果
x=span(1):h:span(2);
n=length(x);
y=zeros(1,n);%存放数值的解
x(1)=x0;
y(1)=y0;
for i=1:n-1
    y(i+1)=y(i)+h.*feval(dufun,x(i),y(i));%Euler格式
end

２．改进的Euler格式

function[x,y]=Gaijin_Euler(func,span,h,x0,y0)
%改进的Euler格式；
%span为区间
%y(x0)=y0;
%n表示区间等份；
x=span(1):h:span(2);
n=length(x);
y=zeros(1,n);
x(1)=x0;
y(1)=y0;
for i=1:n
    y1=y(i)+h*feval(func,x(i),y(i));
    y2=y(i)+h*feval(func,x(i+1),y1);
    y(i)=(y1+y2)/2;
end

 方法二：常见的Rungerkutta3与Rungerkutta4格式

 

１.Rungerkutta3
function [x,y]=Rungekutta3(f1,span,h,x0,y0)
x=span(1):h:span(2);
n=length(x);
y=zeros(1,n);
x(1)=x0;
y(1)=y0;
for j=1:n-1
K1=feval(f1,x(j),y(j));
K2=feval(f1,x(j)+h/2,y(j)+h/2*K1);
K3=feval(f1,x(j)+h,y(j)-h*K1+h*2*K2);
y(j+1)=y(j)+(h/6)*(K1+4*K2+K3);
end

 

１.Rungerkutta4
function [x,y]=Rungekutta4(f1,span,h,x0,y0)
x=span(1):h:span(2);
n=length(x);
y=zeros(1,n);
x(1)=x0;
y(1)=y0;
for j=1:n-1
K1=feval(f1,x(j),y(j));
K2=feval(f1,x(j)+h/2,y(j)+h/2*K1);
K3=feval(f1,x(j)+h/2,y(j)+h/2*K2);
K4=feval(f1,x(j)+h,y(j)+h*K3);
y(j+1)=y(j)+(h/6)*(K1+2*K2+2*K3+K4);
end