参考教材《数值分析》李乃成.梅立泉
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;
trunc函数怎么切除小数点后几位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

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。