四阶龙格-库塔法求解常微分方程的初值问题-matlab通用程序

clear

clc

format long

m=input('请输入常微分方程的阶数m=');

a=input('请输入x下限a=');

b=input('请输入x上限b=');

h=input('请输入步长h=');

ym=input('令y(1,1)=y,y(2,1)=y’,y(3,1)=y’’...请输入ym=','s'); %输入的时候必须按照这个形式输入y1=y(1,1);

if m==1 %一阶初值问题单独求解mm=(b-a)/h;

y(1,1)=input('请输入在初值点的函数值f(a)=');

x=a;

y11(1)=y(1,1);

for k1=2:(mm+1)

y1=y(1,1);

K(1,1)=h*(eval(ym)); %计算K1

x=x+h/2;

y(1,1)=y1+K(1,1)/2;

y1=y(1,1);

K(1,2)=h*(eval(ym)); %计算K2

x=x;

y(1,1)=y1+K(1,2)/2-K(1,1)/2;

y1=y(1,1);

K(1,3)=h*(eval(ym)); %计算K3

x=x+h/2;

y(1,1)=y1+K(1,3)-K(1,2)/2;

y1=y(1,1);

K(1,4)=h*(eval(ym)); %计算K4

y11(k1)=y11(k1-1)+(K(1,1)+2*K(1,2)+2*K(1,3)+K(1,4))/6;

y(1,1)=y11(k1);

x=a+(k1-1)*h;

end

y11

else %高阶初值问题

mm=(b-a)/h; %一共要求解mm个数据点

for k2=1:m %读取初值条件

fprintf('请输入%d阶导数的初值f(%d)(a)=\\n',(k2-1),(k2-1));

y(k2,1)=input('=');

end

for k2=1:m

y22(1,k2)=y(k2,1); %先把初值保存在矩阵y22(m,n)中,m

表示第几个所求点,n表示第n阶初值

end

x=a;

for k4=2:(mm+1) %求解mm个数据点的循环for k=1:(m-1) %计算K1,包括每一阶的K1 K(k,1)=h*y(k+1,1); %y(k+1,1)中k+1表示第k+1阶,1表示第一个点;K(k,1)中k表示阶数,1表示K1

end

K(m,1)=h*(eval(ym));

x=x+h/2; %求解K1之前,先重新对x和y赋值

for k3=1:m

y(k3,1)=y(k3,1)+K(k3,1)/2;

end

for k=1:(m-1) %计算K2

K(k,2)=h*y(k+1,1);

end

K(m,2)=h*(eval(ym));

x=x;

for k3=1:m

y(k3,1)=y(k3,1)-K(k3,1)/2+K(k3,2)/2;

end

for k=1:(m-1) %计算K3

K(k,3)=h*y(k+1,1);

end

K(m,3)=h*(eval(ym));

x=x+h/2;

for k3=1:m

y(k3,1)=y(k3,1)+K(k3,3)-K(k3,2)/2; %这里容易出错

end

for k=1:(m-1) %计算K4

K(k,4)=h*y(k+1,1);

end

K(m,4)=h*(eval(ym));

for k5=1:m

y22(k4,k5)=y22(k4-1,k5)+(K(k5,1)+2*K(k5,2)+2*K(k5,3)+K(k5,4))/6;

%这里,除了要求出下一个点的数值,还要求出相应的导数值

end

for k6=1:m %除了对y(1,1)重新赋值外,还要对y(2,1)等重新赋值

y(k6,1)=y22(k4,k6);

end

x=a+(k4-1)*h;

end

y22(:,1) end