博客

%书籍：常用数值算法及其matlab实现

%第10章 常微分方程初值问题的数值解法,例10.14使用

%四阶龙格库塔方法

function [t,z] = rk4symeq(fun, t0, tf, Za, h)

%fun:微分方程的右表达式

%t0, tn为区间

%Za为初值,是列向量

M = floor(tf-t0)/h ;      %离散点的个数M+1

if t0 >= tf

printf(‘左端点必须小于右端点’);

return;

end

N = length(Za);           %获得变量个数,N

z = zeros(M+1, N);

t = zeros(M +1, 1);

t =[t0 : h :tf]’;

z(1,:) = Za’;            %假设Za为列向量，与微分方程中的变量方向统一，变成行向量

for i = 1:M

K1 =  feval(fun, t(i) , z(i,:));                    %K是行向量

K2 =  feval(fun, t(i)+1/2*h ,z(i,:)+1/2* h*K1);

K3 =  feval(fun, t(i)+1/2*h ,z(i,:)+1/2* h*K2);

K4 =  feval(fun, t(i)+ h ,z(i,:)+ h*K3);

z(i+1,:) = z(i,:) +h/6 *(K1 + 2*K2 + 2*K3 + K4);

end

微分方程如下：

%书籍：常用数值算法及其matlab实现

%第10章 常微分方程初值问题的数值解法

%四阶龙格库塔，例10.14

function s = ex10_4(t,z)

%z是个向量,1*2

%输出s也是向量,1*2

s1 = -5* z(1) + 17*z(2);

s2 = -2*z(1) + 5*z(2) +sin(5*t);

s = [s1 s2];          %输出s为列向量

主函数如下：

clear all;clc;close all;

%书籍：常用数值算法及其matlab实现

%第10章 常微分方程组初值问题的数值解法

%四阶龙格库塔，例10.14

%函数原型 function [t,z] = rk4symeq(fun, t0, tn, Za, h)

format long

t0 = 0; tf = 20;

Za = [1 ; -3];     %x0 = 6; y0 = 4

h1 = 0.1;

h2 = 0.02;

[t1,z1] = rk4symeq(@ex10_4, t0, tf , Za, h1);

figure(1)

plot(t1,z1(:,1),’b’,t1,z1(:,2), ‘r*’)

legend(‘y1′,’y2’)

figure(2)

plot(z1(:,1),z1(:,2));

xlabel(‘x’);ylabel(‘y’);

[t2,z2] = rk4symeq(@ex10_4, t0, tf , Za, h2);

figure(3)

plot(t2,z2(:,1),’b’,t2,z2(:,2), ‘r*’)

legend(‘y1′,’y2’)