Numerical Analysis - 2nd Order Taylor Method to Solve IVP Problem With Code

51 views (last 30 days)
function [y1] = TaylorMethod(f,inter,y0,k)
t(1) = inter(1); %initial time
y1(1) = y0; y2(1) = 0; %setting the initial condition
h = 0.1*(2.^-k); % setting the step size
n = 1/h; % number of steps in terms of the step size.
Ft = @(t,y) diff(f,t);
Fy = @(t,y) diff(f,y);
for i = 1:n
t(i+1) = t(i) + h;
euler_y(i+1) = y(i) + h*f(t(i),y(i));
y1(i+1) = euler_y(i+1) + ((h^2)/2)*(Ft(t(i),y1(i)) + Fy(t(i),y1(i))*f(t(i),y1(i)));
end
Hello,
I am just playing around with some numerical methods and I seem to be having some issues with my code; I want my code to be as robust as possible, so I allocated the variables Fy and Ft to be the partial derivatives of a function 'f' inputed in the function TaylorMethod. However, I get an erro trying to do this.
The IVP I am trying to solve with my function is
y' = 5*(t^4)*y --> this is function 'f'
inter = [0 1]
y(0) = 1
Thank you in advance
  4 Comments
Maria Raheb
Maria Raheb on 10 Apr 2020
@darova
Thank you for the correction, but the problem is arising in my use of function handles for the partial derivatives of my function f.

Sign in to comment.

Answers (1)

Guru Mohanty
Guru Mohanty on 14 Apr 2020
I understand, you want to solve the IVP with 2ndorder Tayler Method. In your code the error occurs during computation of partial differentiation. You can also do the Partial differentiation using Symbolic toolbox and evaluate its value using subs function. Here is a sample code for it.
clc; clear all;
f =@(t,y)(5*(t^4)*y);% --> this is function 'f'
inter = [0 1];
y0 = 1;
k=1;
Out = TaylorMethod(f,inter,y0,k);
plot(Out);
function [y1] = TaylorMethod(f,inter,y0,k)
syms t y tp yp
Ft1=subs(diff(f,t),{t,y},{tp,yp});
Fy1=subs(diff(f,y),{t,y},{tp,yp});
clear t y;
t(1) = inter(1); %initial time
y1(1) = y0;
y2(1) = 0; %setting the initial condition
h = 0.1*(2.^-k); % setting the step size
n = 1/h; % number of steps in terms of the step size.
for i = 1:n-1
t(i+1) = t(i) + h;
euler_y(i+1) = y1(i) + h*f(t(i),y1(i));
y1(i+1) = euler_y(i+1) + ((h^2)/2)*(subs(Ft1,{tp,yp}, {t(i),y1(i)}) + subs(Fy1,{tp,yp},{t(i),y1(i)})*f(t(i),y1(i)));
end
end

Products


Release

R2018b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!