6 MATLAB解微分方程問題 龍格庫塔法

2022-11-24 21:46:49 字數 3083 閱讀 3616

老師說系統給的ode45好多都解決不了。

1.lorenz系統

test.m

clear; 

clc;

%系統龍格庫塔法

[t,h] = ode45(@test_fun,[0

40],[1240

]);plot3(h(:,

1),h(:,2),h(:,3

));grid on;

%自定義龍格庫塔法

[t1,h1]=runge_kutta(@test_fun,[12

40],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];

執行結果(老師給的引數太亂了,調不出來目標的那種樣子):