老师说系统给的ode45好多都解决不了。
1.lorenz系统
test.m
clear; clc; %系统龙格库塔法 [t,h] = ode45(@test_fun,[0 40],[12 4 0]); plot3(h(:,1),h(:,2),h(:,3)); grid on; %自定义龙格库塔法 [t1,h1]=runge_kutta(@test_fun,[12 4 0],0.01,0,40); figure; plot3(h1(1,:),h1(2,:),h1(3,:),'r') grid on;
runge_kutta.m
%参数表顺序依次是微分方程组的函数名称,初始值向量,步长,时间起点,时间终点(参数形式参考了ode45函数) function [x,y]=runge_kutta(ufunc,y0,h,a,b) n=floor((b-a)/h); %步数 x(1)=a; %时间起点 y(:,1)=y0; %赋初值,可以是向量,但是要注意维数 for i=1:n %龙格库塔方法进行数值求解 x(i+1)=x(i)+h; k1=ufunc(x(i),y(:,i)); k2=ufunc(x(i)+h/2,y(:,i)+h*k1/2); k3=ufunc(x(i)+h/2,y(:,i)+h*k2/2); k4=ufunc(x(i)+h,y(:,i)+h*k3); y(:,i+1)=y(:,i)+h*(k1+2*k2+2*k3+k4)/6; end
test_fun.m
%构造微分方程 function dy=test_fun(t,y) a = 16; b = 4; c = 45; dy=[a*(y(2)-y(1)); c*y(1)-y(1)*y(3)-y(2); y(1)*y(2)-b*y(3)];
2.作业一
这里取n = 2
代码:
main.m
clear; clc; %自定义龙格库塔法 %P = [0.1, 0.25, 0.25, 0.25, 0.25]; P = [1.9, 0.8, 0.8, 0.8, 0.8]; [t ,h ]=runge_kutta(@test_fun,P,0.1,0,72); plot(t,h(1,:)); hold on; plot(t,h(1,:)); plot(t,h(2,:)); plot(t,h(3,:)); plot(t,h(4,:)); plot(t,h(5,:)); hold off;
runge_kutta.m
%参数表顺序依次是微分方程组的函数名称,初始值向量,步长,时间起点,时间终点(参数形式参考了ode45函数) function [x,y]=runge_kutta(ufunc,y0,h,a,b) n=floor((b-a)/h); %步数 %x = zeros(1,n) %时间起点 x(1) = a; y(:,1)=y0; %赋初值,可以是向量,但是要注意维数 for i=1:n %龙格库塔方法进行数值求解 x(i+1) = x(i)+h; k1 = ufunc(x(i),y(:,i)); k2 = ufunc(x(i)+h/2,y(:,i)+h*k1/2); k3 = ufunc(x(i)+h/2,y(:,i)+h*k2/2); k4 = ufunc(x(i)+h,y(:,i)+h*k3); y(:,i+1) = y(:,i)+h*(k1+2*k2+2*k3+k4)/6; end
test_fun.m
%构造微分方程 function dy=test_fun(t,P) vs = 0.76, vm = 0.65, vd = 0.95; k1 = 1.9, k2 = 1.3 ,ks = 0.38; M = P(1), P0 = P(2), P1 = P(3), P2 = P(4), PN = P(5); V1 =3.2 , V2 = 1.58, V3 = 5 , V4 = 2.5; K1 = 1, K2 = 2, K3 = 2, K4 = 2, Kd = 0.2, Km = 0.5; n = ; dy= [vs*(K1^n/(K1^n+PN^n))-vm*M/(Km+M); ks*M-V1*P0/(K1+P0)+V2*P1/(K2+P1); V1*P0/(K1+P0)-V2*P1/(K2+P1)-V3*P1/(K3+P1)+V4*P2/(K4+P2); V3*P1/(K3+P1)-V4*P2/(K4+P2)-k1*P2+k2*PN-vd*P2/(Kd+P2); k1*P2-k2*PN];
运行结果(老师给的参数太乱了,调不出来目标的那种样子):