## Help with heat equation

Asked by Janvier Solaris

### Janvier Solaris (view profile)

on 14 Jun 2018
Latest activity Commented on by Janvier Solaris

### Janvier Solaris (view profile)

on 16 Jun 2018
Hello I have the following question I am trying to understand how to do from the heat equation part Can anyone help?

### Tags

Answer by Image Analyst

### Image Analyst (view profile)

on 15 Jun 2018

Try
sigma = 0.2;
r = 0.05;
S0 = 10 : 10 : 100;
K = 50
T = 2
xmax = 10
xmin = -10
N = [10, 100, 200, 400, 800];
M = N
for k = 1 : length(N)
thisM = M(k);
thisN = N(k);
blsprice(................
absoluteError = .........
plot(S0, absoluteError, 'b*-');
hold on;
grid on;
end
I don't know what "the Call Option" or the heat equation or blsprice() is so I can't help with that, but at least this is a start. Presumably blsprice() uses those parameters as inputs. I'm also not sure how to calculate the absolute error.

Janvier Solaris

### Janvier Solaris (view profile)

on 15 Jun 2018
This is what I have so far:
S0= 10:10:100;
K = 50;
r = .05;
T = 2;
sigma = .2;
callV0 = blsprice(S0, K, r, T, sigma)
M = [ 10 100 200 400 800];
N = M;
xmin = -10;
xmax = 10;
dx = (xmax-xmin)/(M+1);
dt = T/N;
r= dt/(dx*dx);
v = zeros(M,N);
x = linspace(xmin,xmax, M)
v(:,1) = max(S0*exp((r-sigma**2/2)*T+sigma*x)-K,0)
%v(:,1) = max(exp(x)-1,0);
for n = 2:N
%lines below are initial and boundary conditions
v(1,n) = v(1, 1);
v(M,n) = v(M,1);
v(2:M-1,n) = (1-2*r)*v(2:M-1,n-1)+r*v(1:M-2,n-1)+r*v(3:M,n-1);
end
error = abs(callV0 - V)
surf(v)
The issue and where I am getting stuck are my boundary conditions. I am assuming xmin and xmax are the boundary conditions but I don't know how to incorporate them where I have initial conditions
Image Analyst

### Image Analyst (view profile)

on 15 Jun 2018
You can't do
for n = 2:N
because N is an array. You have to think what you really want. You need to extract the values of M and N from the array and use them somehow.
Even if you did
for n = N
which is allowed, n would take on the value of N(1) first, which is 10. So then you say
v(1, 10) = v(1,1);
well, what about the rest of the array? What about columns 2 through 9?
You also shouldn't do
v = zeros(M,N);
because they're arrays. You can if you use thisM and thisN and put that inside the loop as I showed you. Take another crack at it.
Janvier Solaris

### Janvier Solaris (view profile)

on 16 Jun 2018
Hi @Image Analyst, I was able to modify my code as such which works similarly to your modifications;
S0 = 10:10:100;
K = 50;
r = .05;
T = 2;
sigma = .2;
price = blsprice(S0,K,r,T,sigma)
amax = 10;
amin = -10;
M = [10 100 200 400 800];
N = [10 100 200 400 800];
for k = 1:5
dx = (amax - amin)/(M(k)-1);
dt = T/(N(k)-1);
v = zeros(M(k), N(k));
r= dt/dx^2;
x = linspace(amin,amax,M(k));
v(:,1) = max(exp(x)-1,0);
for n = 2:N(k)
%2 lines below are initial and boundary conditions
v(1,n) = v(1, 1);
v(N(k),n) = v(M(k),1);
v(2:M(k)-1,n) = (1-2*r)*v(2:M(k)-1,n-1)+r*v(1:M(k)-2,n-1)+r*v(3:M(k),n-1);
end
approx = v(2:M(k)-1,n)
end
error = abs(approx - price)
B= table(approx')
G=table(price')
surf(v)
the loop grabs every item within the array as thisN and thisM would have .
Now the only issue I am getting is that I am getting the error, the limits are too large and it is not plotting for my 3d plot. when I change the array to for k = 1: 3 then it works but not for 400 and higher would you know anything about "limits are too large" error?