微分方程數值解法的matlab程式

微分方程數值解法的matlab程式

尤拉法

clear
clc
syms x y
f=2*x;
a=0;
b=10;
y0=0;
h=0.1;
n=(b-a)/h+1;
xx=zeros(n,1);
yy=zeros(n,1);
for i=1:n
    xx(i)=a+h*(i-1); 
end
yy=zeros(n,1);
yy(1)=y0;
for i=2:n
    yy(i)=yy(i-1)+h*subs(f,{x,y},{xx(i-1),yy(i-1)}); 
end
ff1=[xx,yy];

改進尤拉法

clear
clc
syms x y
f=3*x^2;
a=0;
b=10;
y0=0;
h=0.1;
n=(b-a)/h+1;
xx=zeros(n,1);
yy=zeros(n,1);
for i=1:n
    xx(i)=a+h*(i-1); 
end
yy=zeros(n,1);
yy(1)=y0;
for i=2:n
    yy(i)=yy(i-1)+h*subs(f,{x,y},{xx(i-1),yy(i-1)});
    p=zeros(4,1);
    p(1)=yy(i);
    for j=2:4
       p(j)=yy(i-1)+(h/2)*(subs(f,{x,y},{xx(i-1),yy(i-1)})+subs(f,{x,y},{xx(i),p(j-1)}));
    end
    yy(i)=p(4);
end
ff1=[xx,yy];

龍格庫塔法

clear
clc
syms x y
f=3*x^2;
a=0;
b=10;
y0=0;
h=0.1;
n=(b-a)/h+1;
xx=zeros(n,1);
yy=zeros(n,1);
for i=1:n
    xx(i)=a+h*(i-1); 
end
yy=zeros(n,1);
yy(1)=y0;
for i=2:n
    k1=subs(f,{x,y},{xx(i-1),yy(i-1)});
    k2=subs(f,{x,y},{xx(i-1)+h/2,yy(i-1)+h/2*k1});
    k3=subs(f,{x,y},{xx(i-1)+h,yy(i-1)-h*k1+2*h*k2});
    yy(i)=yy(i-1)+h/6*(k1+4*k2+k3);
end
ff1=[xx,yy];

隱式龍格庫塔法

%隱式比顯式有更高的精度，但是因為求解非線性或者線性方程組導致計算量過大。
%下面編寫三階隱式
clear
clc
syms x y
f=3*x^2;
a=0;
b=10;
y0=0;
h=0.1;
n=(b-a)/h+1;
xx=zeros(n,1);
yy=zeros(n,1);
for i=1:n
    xx(i)=a+h*(i-1); 
end
yy=zeros(n,1);
yy(1)=y0;
c=[(5-15^0.5)/10,0.5,(5+15^0.5)/10];
b=[5/18,4/9,5/18];
aa=[5/36,(10-3*15^0.5)/45,(25-6*15^0.5)/180;(10+3*15^0.5)/72,2/9,(10-3*15^0.5)/72;(25+6*15^0.5)/180,(10+3*15^0.5)/45,5/36];
for i=2:n
    syms k1 k2 k3
    k1=subs(f,{x,y},{xx(i-1)+c(1)*h,yy(i-1)+h*(aa(1,1)*k1+aa(1,2)*k2+aa(1,3)*k3)});
    k2=subs(f,{x,y},{xx(i-1)+c(2)*h,yy(i-1)+h*(aa(2,1)*k1+aa(2,2)*k2+aa(2,3)*k3)});
    k3=subs(f,{x,y},{xx(i-1)+c(3)*h,yy(i-1)+h*(aa(3,1)*k1+aa(3,2)*k2+aa(3,3)*k3)});
    yy(i)=yy(i-1)+h/6*(k1+4*k2+k3);
end
ff1=[xx,yy];

半隱式龍格庫塔法

%半隱式比隱式計算量更小，並且精度相當。
%下面編寫二階半隱式
clear
clc
syms x y
f=6*x;
A=diff(f,y,1);
a=0;
b=10;
y0=0;
h=0.1;
n=(b-a)/h+1;
xx=zeros(n,1);
yy=zeros(n,1);
for i=1:n
    xx(i)=a+h*(i-1); 
end
yy=zeros(n,1);
yy(1)=y0;
c=[(5-15^0.5)/10,0.5,(5+15^0.5)/10];
b=[5/18,4/9,5/18];
aa=[5/36,(10-3*15^0.5)/45,(25-6*15^0.5)/180;(10+3*15^0.5)/72,2/9,(10-3*15^0.5)/72;(25+6*15^0.5)/180,(10+3*15^0.5)/45,5/36];
omig1=(6+6^0.5)/6;
omig2=(6-6^0.5)/6;
c2=0.17378667;
b1=-0.41315432;
b2=1-b1;
for i=2:n
    k1=(1-omig1*subs(A,{x,y},{xx(i-1),yy(i-1)}))^(-1)*subs(f,{x,y},{xx(i-1),yy(i-1)});
    k2=(1-omig2*subs(A,{x,y},{xx(i-1)+c2*h,yy(i-1)+h*c2*k1}))^(-1)*subs(f,{x,y},{xx(i-1)+c2*h,yy(i-1)+h*c2*k1});
    yy(i)=yy(i-1)+h*(b1*k1+b2*k2);
end
ff1=[xx,yy];