Math 563 Scientific Computing II

Homework 2 Due Jan.17 Hangjie Ji

Problem:

Write a program to interpolate the function f(x,y) = cos(x) sin(y) at 16 equally-spaced points in the square (-pi,pi) x (-pi,pi), using multi-dimensional polynomials that are cubics in each coordinate. Plot the error in the function value, and discuss how the error compares with the theory in section 10.3.4.

Solution: Generalizing the one dimensional Lagrange interpolation, we obtain the following method. Given 16 equally-spaced points (x_i,y_j) for i,j=0,1,2,3 and the corresponding function values f(x_i,y_j), for (x,y)\in (-\pi,\pi)\times(-\pi,\pi),the basis polynomials are

\phi_{i,j}(x,y)=\prod_{m\neq i}\frac{x-x_m}{x_i-x_m}\prod_{n\neq j}\frac{y-y_n}{y_j-y_n}.

In this case, the coefficients of these basis polynomials in the interpolaton are f(x_i,y_j). Therefore the resulting Lagrange interpolation is

f(x,y)\approx p(x,y)=\sum_{i=0}^3\sum_{j=0}^3 f(x_i,y_j)\phi_{i,j}(x,y).

The interpolation result is shown in Figure 1, where the red dots represent the interpolation points. Here I evaluated the function values at 80*80 equally spaced points in the domain. The error in the function value is plotted in Figure 2, with the maximum error reaching 0.4359.

According to Theorem 10.3.2 in section 10.3.4, we have

\sup_{(x,y)\in \Omega}|f(x,y)-p(x,y)|\le h^4(\sum_{|\alpha|=4}\frac{1}{\alpha!}\sup_{(x,y)\in \Omega}|D^{\alpha}f(x,y)|)(\sum_{i=0}^{3}\sum_{j=0}^3\sup_{(x,y)\in \Omega}|\phi_{i,j}(x,y)|}).

Since f(x,y)=cos(x)sin(y), we have \sup_{(x,y)\in \Omega}|D^{\alpha}f(x,y)|\le 1,\forall \alpha.

Since the domain has diameter h=2pi, and \sup_{(x,y)\in \Omega}|\phi_{i,j}(x,y)|=1,\forall i,j=0,1,2,3.,

We conclude that

\sup_{(x,y)\in \Omega}|f(x,y)-p(x,y)|\le (2\pi)^4 \times \sum_{|\alpha|=4}\frac{1}{\alpha !}\times 16 = \frac{704}{3}\pi^4.

This is not a good error estimate in this case, but at least the error in numerical result satisfies the condition.

Figure 1

Figure 2

The Matlab code is attached here: hw2
% interpolate the function f(x,y) = cos(x) sin(y) at 16 equally-spaced
% points in the square (-pi,pi) x (-pi,pi), using multi-dimensional polynomials
% that are cubics in each coordinate

close all;
clear all;
clc;

n = 80;
[x_ip,y_ip] = meshgrid(-pi:pi*2/3:pi,-pi:pi*2/3:pi);
f = cos(x_ip).*sin(y_ip);
x_ip = linspace(-pi,pi,4); % interpolation points
y_ip = linspace(-pi,pi,4);
f = zeros(4,4);
for i = 1 : 4
    for j = 1 : 4
        f(i,j) = cos(x_ip(i)).* sin(y_ip(j));
    end
end

x = linspace(-pi,pi,n);
y = linspace(-pi,pi,n);
p = zeros(n,n); % interpolation function
error = zeros(n,n); % error
f_true = zeros(n,n); % exact function;

for i = 1:n
    for j = 1:n
        f_estimate = 0;
        for k = 1 : 4
            for s = 1 : 4
            temp1 = [x_ip(1:k-1),x_ip(k+1:4)];
            temp2 = [y_ip(1:s-1),y_ip(s+1:4)];
            f_estimate = f_estimate + f(k,s)*prod(x(i)-temp1)/prod(x_ip(k)-temp1)...
                *prod(y(j)-temp2)/prod(y_ip(s)-temp2);
            end
        end
        p(i,j) = f_estimate;
        f_true(i,j) = cos(x(i))*sin(y(j));
        error(i,j) = f_true(i,j) - p(i,j);
    end
end


figure(1),mesh(x,y,p'),xlabel('x','FontSize',14),ylabel('y','FontSize',14);
hold on;
for i = 1:4
    for j = 1:4
    plot3(x_ip(i),y_ip(j),f(i,j),'ro');
    end
end
hold off;
figure(2),mesh(x,y,error'),xlabel('x','FontSize',14),ylabel('y','FontSize',14),
title('error','FontSize',14);
err_max = max(max(error));


Published with MATLAB® 7.9

Last modified: